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Abstract  Spin-orbit  coupling  can  be  described  in  two  approaches.  The  first  method,  known 
as  the  “MacDonald  torque  ”,  is  often  combined  with  a  convenient  assumption  that  the  quality 
factor  Q  is  frequency-independent.  This  makes  the  method  inconsistent,  because  derivation 
of  the  expression  for  the  MacDonald  torque  tacitly  fixes  the  rheology  of  the  mantle  by  making 
Q  scale  as  the  inverse  tidal  frequency.  Spin-orbit  coupling  can  be  treated  also  in  an  approach 
called  “ the  Darwin  torque  While  this  theory  is  general  enough  to  accommodate  an  arbitrary 
frequency-dependence  of  Q,  this  advantage  has  not  yet  been  fully  exploited  in  the  literature, 
where  Q  is  often  assumed  constant  or  is  set  to  scale  as  inverse  tidal  frequency,  the  latter  asser¬ 
tion  making  the  Darwin  torque  equivalent  to  a  corrected  version  of  the  MacDonald  torque. 
However  neither  a  constant  nor  an  inverse-frequency  Q  reflect  the  properties  of  realistic 
mantles  and  crusts,  because  the  actual  frequency-dependence  is  more  complex.  Hence  it  is 
necessary  to  enrich  the  theory  of  spin-orbit  interaction  with  the  right  frequency-dependence. 
We  accomplish  this  programme  for  the  Darwin-torque-based  model  near  resonances.  We 
derive  the  frequency-dependence  of  the  tidal  torque  from  the  first  principles  of  solid-state 
mechanics,  i.e.,  from  the  expression  for  the  mantle’s  compliance  in  the  time  domain.  We 
also  explain  that  the  tidal  torque  includes  not  only  the  customary,  secular  part,  but  also  an 
oscillating  part.  We  demonstrate  that  the  Impq  term  of  the  Darwin-Kaula  expansion  for  the 
tidal  torque  smoothly  passes  zero,  when  the  secondary  traverses  the  Impq  resonance  (e.g., 
the  principal  tidal  torque  smoothly  goes  through  nil  as  the  secondary  crosses  the  synchronous 
orbit).  Thus,  we  prepare  a  foundation  for  modeling  entrapment  of  a  despinning  primary  into 
a  resonance  with  its  secondary.  The  roles  of  the  primary  and  secondary  may  be  played,  e.g., 
by  Mercury  and  the  Sun,  correspondingly,  or  by  an  icy  moon  and  a  Jovian  planet.  We  also 
offer  a  possible  explanation  for  the  “improper”  frequency-dependence  of  the  tidal  dissipation 
rate  in  the  Moon,  discovered  by  LLR. 
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1  Introduction 

The  paper  addresses  tidal  torques  acting  on  a  body  in  the  vicinity  of  a  spin-orbit  resonance. 
Although  the  topic  has  been  popular  since  mid-sixties  and  has  already  been  addressed  in 
books,  the  common  models  are  not  entirely  adequate  to  the  actual  physics.  Just  as  in  the 
nonresonant  case  discussed  in  Efroimsky  and  Williams  (2009),  a  generic  problem  with  the 
popular  models  of  libration  or  of  capture  into  a  resonance  is  that  they  employ  wrong  rhe¬ 
ologies  (the  work  by  Rambaux  et  al.  2010  being  the  only  exception  we  know  of).  Above 
that,  the  model  based  on  the  MacDonald  torque  suffers  a  defect  stemming  from  a  genuine 
inconsistency  inherent  in  the  theory  by  MacDonald  (1964)  and  Gerstenkorn  (1955). 

As  explained  in  Efroimsky  and  Williams  (2009)  and  Williams  and  Efroimsky  (2012),  the 
MacDonald  theory,  both  in  its  original  and  corrected  versions,  tacitly  fixes  an  unphysical 
shape  of  the  functional  dependence  Q(x),  where  Q  is  the  dissipation  quality  factor  and  /  is 
the  tidal  frequency.  So  we  base  our  treatment  on  the  developments  by  Darwin  (1879,  1880) 
and  Kaula  (1964),  combining  those  with  a  realistic  frequency-dependence  of  the  damping 
rate. 

Since  our  main  purpose  is  to  lay  the  groundwork  for  the  subsequent  study  of  the  process 
of  falling  into  a  resonance,  the  two  principal  results  obtained  in  this  paper  are  the  following: 

(a)  Starting  with  a  realistic  rheological  model  (an  expression  for  the  compliance  in  the  time 
domain),  we  derive  the  complex  Love  numbers  kt  as  functions  of  the  frequency  x,  and 
write  down  their  negative  imaginary  parts  as  functions  of  the  frequency:  —2m  [kj  (x )]  = 
\kj  (x ) I  sin  (x).  It  is  these  expressions  that  appear  as  factors  in  the  terms  of  the  Dar- 
win-Kaula  expansion  of  tides.  These  factors’  frequency-dependencies  demonstrate  a 
nontrivial  shape,  especially  near  resonances.  This  shape  plays  a  crucial  role  in  mod¬ 
eling  of  despinning  in  general,  specifically  in  modeling  the  process  of  falling  into  a 
spin-orbit  resonance. 

(b)  We  demonstrate  that,  beside  the  customary  secular  part,  the  Darwin  torque  contains  a 
usually  omitted  oscillating  part. 

A  much  extended  version  of  this  paper  is  available  online  (arXiv:  1 105.6086). 


2  Linear  bodily  tides 

Linearity  of  tide  means  that:  (a)  under  a  static  load,  deformation  scales  linearly,  and  (b)  under 
undulatory  loading,  the  same  linear  law  applies,  separately,  to  each  frequency  mode.  Thus 
the  deformation  magnitude  at  each  frequency  should  depend  linearly  upon  the  tidal  stress 
at  this  frequency,  and  should  bear  no  dependence  upon  loading  at  other  tidal  modes.  So  the 
dissipation  rate  at  that  frequency  will  depend  on  the  stress  at  that  frequency  only. 

2. 1  Linearity  of  the  tidal  deformation 

At  a  point  R  =  (R,X,<j>),  the  potential  due  to  a  tide-raising  secondary  of  mass  M*ec, 
located  at  r*  =  (r*,  A*,  (j>*)  with  r*  >  R,  is  expandable  over  the  Legendre  polynomials 
Pi  (cosy): 
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W(R,  r*)  =  X  W,( R,  r*)  =  X  (£  ) 

1=2  1=2  V  7 

GM*ec  /  R  \ 

- — ^Z(^)  z 

1=2  X  7  m= 0 


—  )  P/(cosy) 


(l  —  m)\ 
(/  +  m)\ 


(2  -  <50m  )/’/,„  (sin  4>)Pim(s'm(f>*) 


x  cos/ra(A  —  A*), 


(1) 


where  G  =  6.7  x  10_n  m3kg^1  s~2  is  Newton’s  gravity  constant,  and  y  is  the  angular 
separation  between  the  vectors  r*  and  R  pointing  from  the  primary’s  centre.  The  latitudes 
4>,  cj>*  are  reckoned  from  the  primary’s  equator,  while  the  longitudes  A ,  A*  are  reckoned  from 
a  fixed  meridian. 

Under  the  assumption  of  linearity,  the  term  W/(R,  r*)  in  the  secondary’s  potential  causes 
a  linear  deformation  of  the  primary’s  shape.  The  subsequent  adjustment  of  the  primary’s 
potential  being  linear  in  the  said  deformation,  the  adjustment  U/  of  this  potential  is  propor¬ 
tional  to  W] .  The  theory  of  potential  requires  U\  (r)  to  fall  off,  outside  the  primary,  as  +1^ . 
Thus  the  overall  amendment  to  the  potential  of  the  primary  amounts  to: 


°°  °°  /  R\l+l 

U (r)  =  X  Ut  (r)  =  X  k‘  (  -  )  W,  (R,  r*) 
1=2  1=2  '  r  ' 


(2) 


R  now  being  the  mean  equatorial  radius  of  the  primary,  R  =  ( R .  </>,  A)  being  a  surface  point, 
r  =  (r,  <p,  A)  being  an  exterior  point  located  above  it  at  a  radius  r  >  R.  The  coefficients  k/, 
called  Love  numbers,  are  defined  by  the  primary’s  rheology. 

For  a  homogeneous  incompressible  spherical  primary  of  density  p,  surface  gravity  g,  and 
rigidity  /x,  the  static  Love  number  of  degree  /  is  given  by 


k, 


3  1 

2(1  -  1)  1  +  A/’ 


where  A / 


(2  /“  — F  4/  — L  3 )/x 

IgpR 


3  (2  Z2  +  4  /  +  3)  p, 
4/jrGp2  R2 


For  R  r,r*,  consideration  of  the  1  =  2  input  in  (2)  turns  out  to  be  sufficient.1 

These  formulae  apply  to  static  deformations.  Flowever  an  actual  tide  is  never  static,  except 
in  the  case  of  synchronous  orbiting  with  a  zero  eccentricity  and  inclination.2  A  realistic  per¬ 
turbing  potential  produced  by  the  secondary  carries  a  spectrum  of  modes  a>lm  (positive 
or  negative)  numbered  with  four  integers  linpq  as  in  formula  (101)  below.  The  perturbation 
causes  a  spectrum  of  stresses  in  the  primary,  at  frequencies  Ximpq  =  \wimpq  I-  Although  in 
a  linear  medium  strains  are  generated  at  the  frequencies  of  the  stresses,  friction  makes  each 
Fourier  component  of  the  strain  fall  behind  the  appropriate  component  of  the  stress.  Friction 
also  reduces  the  magnitude  of  the  shape  response — hence  the  deviation  of  a  dynamical  Love 
number  A/(x)  from  its  static  counterpart  ki  =  k/( 0).  Below  we  shall  explain  that  formulae 
(2-3)  can  be  easily  adjusted  to  the  case  of  undulatory  tidal  loads  in  a  homogeneous  planet 
or  in  tidally-despinning  homogeneous  satellite  (treated  now  as  the  primary,  with  its  planet 
playing  the  role  of  the  tide-raising  secondary).  Flowever  generalisation  of  formulae  (2-3)  to 
the  case  of  a  librating  moon  (treated  as  a  primary)  turns  out  to  be  nontrivial,  because  the 


1  Special  is  the  case  of  Phobos,  for  whose  orbital  evolution  the  k 3  and  perhaps  even  the  F4  terms  may  be 
relevant  (Bills  et  al.  2005).  Another  class  of  exceptions  is  constituted  by  close  binary  asteroids.  The  topic  is 
addressed  by  Taylor  and  Margot  (2010),  who  took  into  account  the  Love  numbers  up  to  k 

2  The  case  of  a  permanently  deformed  moon  in  a  1:1  spin-orbit  resonance  falls  under  this  description  too. 
Recall  that  in  the  tidal  context  the  distorted  body  is  taken  to  be  the  primary.  So  from  the  viewpoint  of  the 
satellite  its  host  planet  is  orbiting  the  satellite  synchronously,  thus  creating  a  static  tide. 
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standard  derivation  by  Love  (1909,  1911)  falls  apart  in  the  presence  of  the  non-potential 
inertial  force  containing  the  time-derivative  of  the  primary’s  angular  velocity. 

The  frequency-dependence  of  a  dynamical  Love  number  takes  its  origins  in  the  “inertia” 
of  strain  and,  therefore,  of  the  shape  of  the  body.  Hence  the  analogy  to  linear  circuits:  the 
Ith  components  of  W  and  U  act  as  a  current  and  voltage,  while  the  lrh  Love  number  plays, 
up  to  a  factor,  the  role  of  impedance.  Therefore,  under  a  sinusoidal  load  of  frequency  /,  it  is 
convenient  to  replace  the  actual  Love  number  with  its  complex  counterpart 

k(x)  =  \k(x)\ exp  [-/€;(/)],  (4) 

e/  being  the  phase  delay  of  the  reaction  relative  to  the  load  (Munk  and  MacDonald  1960; 
Zschau  1978).  The  “minus”  sign  in  (4)  makes  U  lag  behind  W  for  a  positive  e/.  (So  the 
situation  resembles  a  circuit  with  a  capacitor,  where  the  current  leads  the  voltage.) 

For  a  steady  deformation,  the  lag  should  vanish,  and  so  should  the  entire  imaginary  part: 

Tm  [k,  (0)]  =  |  k,  (0) |  sin  e, (0)  =  0,  (5) 

leaving  the  complex  Love  number  real: 

k,  (0)  =  Tie  [k,  (0)]  =  |  k,  (0)  |  cos  6,  (0) ,  (6) 

and  equal  to  the  customary  static  Love  number: 

ki(0)  =  k,.  (7) 

Combined  with  the  constitutive  (rheological)  law,  the  equation  of  motion  renders  the 
complex  k/(x),  as  explained  in  Appendix  B.l.  Once  k/(x)  is  found,  its  absolute  value 

klix)  =  I*/(X)I  (8) 


and  negative  argument 


6/(x)  =  —  arctan 


1m  [k,(x)\ 
Tie  [fe/(x)] 


(9) 


should  be  inserted  into  the  l,h  term  of  the  Fourier  expansion  for  the  tidal  potential.  Things 
get  simplified  when  we  study  how  the  tide,  caused  on  the  primary  by  a  secondary,  is  acting 
on  that  same  secondary.  In  this  case,  the  lrh  term  in  the  Fourier  expansion  contains  |k/(x)| 
and  6/ (x )  in  the  convenient  combination  k/(x)  sine/(x),  which  is  exactly  —  Tm  [k/(x)]- 
Rigorously  speaking,  we  should  say  not  “the  l,h  term”,  but  “the  ITh  terms”,  as  each  /  cor¬ 
responds  to  an  infinite  set  of  positive  and  negative  Fourier  modes  CL>impq ,  the  physical  forcing 
frequencies  being  x  =  Xhnpq  =  \<*>lmPq\-  Thus,  while  the  functional  forms  of  both  |k/(x)| 
and  sin  e/(x)  depend  only  on  /,  these  functions  take  values  that  are  different  for  different  sets 
of  numbers  mpq,  because  x  assumes  different  values  ximpq  on  these  sets.  For  triaxial  bodies 
the  functional  forms  of  |k/(x)|  and  sine/(x)  may  depend  also  on  m,  p,  q. 


2.2  Damping  of  a  linear  tide 

Besides  the  standard  assumption  f//(r)  oc  W/(R,  r*),  the  linearity  condition  includes  the 
requirement  that  the  functions  k/(x)  and  e,  (x)  be  well  defined.  This  implies  that  they  depend 
solely  upon  the  frequency  x ,  not  upon  the  other  frequencies  involved.  Nor  shall  the  Love 
numbers  or  lags  be  influenced  by  the  stress  or  strain  magnitudes  at  this  or  other  frequencies. 
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Then,  at  frequency  x ,  the  mean  (over  a  period)  damping  rate  (E  (x ))  depends  on  the  value 
of  x  and  on  the  loading  at  that  frequency,  and  is  not  influenced  by  the  other  frequencies: 

<£(X)>  =-XEpeak(x)/Q(x)  (10) 

or,  equivalently: 


A£cyc/e(x)  =  -2?r  Epeak(X)/Q(x),  (11) 


A Ecycie(x)  being  the  one-cycle  energy  loss,  and  Q(x)  being  the  so-called  quality  factor. 

If  E peak  (x )  in  (10-11)  is  agreed  to  denote  the  peak  energy  stored  at  frequency  x>  the 
appropriate  Q  factor  is  connected  to  the  phase  lag  e(x)  through 


Q  =  sin  |e| 

energy 


(12) 


and  not  through  Gen*  =  tan  |e|  as  often  presumed  (see  Appendix  A  for  explanation). 

If  Epeai :(x)  is  defined  as  the  peak  work ,  the  appropriate  Q  factor  is  related  to  the  lag  via 


Q  1  ,  =  tan 

work 


kl  [l  -  (  -  |e|)tan|e|] 


as  demonstrated  in  Appendix  A  below.3  In  the  limit  of  a  small  e,  (13)  becomes 
Qwlk  =sinkl  +  0(e2)  =  kl  +  O(e2), 


(13) 


(14) 


so  definition  (13)  makes  1/  Q  a  good  approximation  to  sin  e  for  small  lags  only. 

For  the  lag  approaching  7r/2,  the  quality  factor  defined  through  (12)  attains  its  minimum, 
Qewrgy  =  *’  whHe  definition  (13)  furnishes  Qwork  =  0.  The  latter  is  not  surprising,  as  in 
the  said  limit  no  work  is  carried  out  on  the  system. 

Linearity  requires  k/(x)  and  therefore  e/(x)  to  be  well-defined  functions,  i.e.,  to  be  inde¬ 
pendent  from  all  the  other  frequencies  but  x-  As  we  see,  the  requirement  extends  to  Q(x). 

The  third  definition,  Q^0idreich  =  tan  kl>  offered  by  Goldreich  (1963),  relates  neither  to 
the  peak  work  nor  to  the  peak  energy.  The  ambiguity  in  definition  of  Q  makes  this  factor 
redundant,  albeit  a  part  of  the  historical  tradition.  Practical  calculations  contain  the  products 
of  the  Love  numbers  by  the  sines  of  the  phase  lags,  k /  sine/,  where  1  is  the  degree  of  the 
appropriate  spherical  harmonic.  One  may  though  define  the  tidal  Q  by  (12)  and  equip  it  with 
the  subscript  /,  thus  emphasising  its  otherness  from  the  seismic  Q  (Efroimsky  2012). 


3  Several  basic  facts  from  continuum  mechanics 


This  section  offers  a  squeezed  synopsis  of  the  basic  linear  solid-state  mechanics.  A  more 
detailed  introduction,  including  a  glossary  and  examples,  is  offered  in  the  Electronic  Sup¬ 
plementary  Material  1  “ Continuum  mechanics.  A  celestial-mechanician’s  survival  kit”. 

3. 1  Stationary  linear  deformation  of  isotropic  incompressible  media 

Mechanical  properties  of  a  medium  are  furnished  by  the  so-called  constitutive  equation  or 
constitutive  law,  which  interrelates  the  stress  tensor  §  with  the  strain  tensor  U  defined  as 

U  =  ^  [(V  ®  u)  +  (V  ®  u/  j  ,  (15) 

u  being  the  displacement.  The  viscous,  elastic  and  hereditary  stresses  are  related  to  U  via 

3  Deriving  this  formula  in  Efroimsky  and  Williams  (2009),  we  inaccurately  termed  Ep eak  (x)  as  peak  energy. 
However  our  calculation  of  Q  was  carried  out  in  understanding  that  Epeak{x )  is  the  peak  work. 
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00  9  «0  TO  . 

§  =  A — U,  §  =  B  U,  §  =  B  U,  (16) 

dt 

where  A  is  a  four-dimensional  matrices  of  viscosities,  B  is  that  of  elastic  moduli,  while  B  is 
an  operator-valued  matrix.  To  furnish  the  value  of  er-  =  ]TW  Bijkl  uki  at  time  A  the  operators 
Bjjki  integrate  the  values  of  ukl  (/)  over  the  interval  t'  G  (—  oo,  t  ]. 

In  isotropic  media,  each  matrix  includes  two  terms  only.  The  elastic  stress  becomes: 

(«>  GO  GO  /I  \  /  1  \ 

§  =  §  volumetric  +  §  devia,oric=  3  K  1  SP  U  j  +  2  p  (U  -  -  I  Sp  Uj  ,  (17) 

with  K  and  fi  being  the  bulk  elastic  modulus  and  the  shear  elastic  modulus,  correspondingly, 
I  standing  for  the  unity  matrix,  and  Sp  denoting  the  trace  of  a  matrix:  Sp  UJ  =  JT  Up . 

The  hereditary  stress  becomes: 

w  to  to  .  n  \  _  /  i  \ 

§  =  §  volumetric  +  §  devia,oric=  3  K  (j  I  Sp  Uj  +  2 /i  (U  -  -  I  Sp  U  j,  (18) 

where  A  and  jl  are  the  bulk-modulus  operator  and  the  shear-modulus  operator,  accordingly. 
The  viscous  stress  acquires  the  form: 

(«)  («)  («)  9/1  \  3/1  \ 

§  =  §  volumetric  +  §  deviatoric=  3  S  ^  !  SP  U  j  +  2  r,  -  \  V  ~  -  I  Sp  U  j  ,  (19) 

the  quantities  f  and  being  termed  as  the  bulk  and  shear  viscosity,  correspondingly. 

The  term  ^  I  Sp  U  is  called  the  volumetric  part  of  the  strain,  while  U  —  -j  I  Sp  U  is  called 
the  deviatoric  part.  Accordingly,  in  expressions  (17-19)  for  the  stresses,  the  pure-trace  terms 
are  called  volumetric,  the  other  term  being  named  deviatoric. 

If  an  isotropic  medium  is  also  incompressible,  the  relative  change  of  the  volume,  Sp  U  =  0, 
vanishes;  so  the  volumetric  parts  of  the  strain  and  stresses  become  nil.  The  incompressibil¬ 
ity  assumption  may  be  applicable  both  to  crusty  and  to  low-porous  icy  bodies.  At  least  for 
Iapetus,  the  low-porosity  assumption  is  likely  to  be  correct  (Castillo-Rogez  et  al.  201 1). 

3.2  Evolving  stresses  and  strains:  basic  notations 

In  the  general  case,  loading  varies  in  time,  so  one  has  to  deal  with  the  stress  and  strain  tensors 
as  functions  of  time.  However,  treatment  of  viscoelasticity  turns  out  to  be  simpler  in  the 
frequency  domain,  i.e.,  in  the  language  of  complex  rigidity  and  complex  compliance.  To  this 
end,  the  stress  cryv  and  strain  uyv  in  a  linear  medium  can  be  Fourier-expanded  as 

OO  OO 

Oyv(t)  =  ^  OyV(Xn)  COs[xnt  +  <pa(xn)]  =  X  ne[dYv(Xn)  exP  (  *XnO]  ■  (2°) 

n=  0  n= 0 

oo  oo 

Uyv(t)  =  X  “yr(x„)  cos[x„t  +  «5„(x„)]  =  X  M“/r(X„)  exP  ('X„0]  -  (21) 

n=  0  n=0 

where  the  complex  amplitudes  are: 

°yv(x)  =  cryv(x)exp['V<r(x)].  i<yv(x)  =  “yv(x)exp['VM(x)]  •  (22) 

while  the  initial  phases  (pa(x)  and  <Pu(x)  are  chosen  in  a  manner  that  sets  the  real  amplitudes 
cryv(x„)  and  uyv(xn)  non-negative. 
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We  wrote  the  above  expansions  as  sums  over  a  discrete  spectrum,  as  the  spectrum  gener¬ 
ated  by  tides  is  discrete.  Generally,  the  sums  can  be  replaced  with  integrals  over  frequency: 

oo  oo 

f  _  i’x*  f  _  ixt 

cryv(t)=  /  cryv(x)e  dx  and  uyv(t)  =  /  uyv(x)  e  dx .  (23) 

0  0 

Whenever  necessary,  the  frequency  is  set  to  approach  the  real  axis  from  below:  Tm(x )  — >  0— . 


3.3  Should  we  consider  positive  frequencies  only? 


At  first  glance,  the  above  question  appears  pointless,  because  in  classical  physics  a  negative 
frequency  is  a  mere  abstraction.  Mathematically,  a  full  Fourier  decomposition  of  a  real  field 
can  always  be  reduced  to  a  decomposition  over  positive  frequencies  only. 

For  example,  the  full  Fourier  integral  for  the  stress  can  be  written  as 


Jyv 


oo  oo 

(t)  =  J  SyV(cj))e,co,dw  =  j  [sYV(x)e'xt  +syv(-x)e~IX‘]  dx,  (24) 

— oo  0 


where  we  define  x  =  \u>\.  Denoting  complex  conjugation  with  asterisk,  we  write: 

-ixt ' 


ayV(t)  = 


0 


s*yv(-x)e,XT  +s*Vv(x)e 


yvv 


dx- 


(25) 


The  stress  is  real:  a*  (f)  =  ayv(t).  Equating  the  right-hand  sides  of  (24)  and  (25),  we  obtain 


■Sxv(-X)  =  Syv(X). 


(26) 


whence 


Vyv(t)  = 


-  ,  ,  ‘X1  -*  ,  .  -ixt 

Syvix)e  +.?  v(x)e 


‘X‘ 


dx  =  Ke  /  2  syv(x)e  dx-  (27) 


This  leads  us  to  (23),  if  we  set 


ffyv(x)  =  2syv(x).  (28) 

Our  switch  from  ayv{t)  =  syv(co)  e  du>  to  the  expansion  cry v(f)  =  J„  OyV(w)e 
dx  makes  things  simpler,  but  the  simplification  comes  at  a  cost,  as  we  shall  see  in  a  second. 
Recall  that  the  tide  can  be  expanded  over  the  modes 

a>impq  =  (/  —  2 p)  6)  +  (l  —  2p  +  q)  M.  +  m  (Si  —  9)  (Z  —  2p  +  q)n  —  m  9,  (29) 


each  of  which  can  assume  positive  or  negative  values,  or  be  zero.  Flere  /,  m,  p,  q  are  some 
integers,  6  is  the  primary’s  sidereal  angle,  9  is  its  spin  rate,  while  u>,  f2,  A4  and  n  are  the 
secondary’s  periapse,  node,  mean  anomaly,  and  mean  motion.  The  tidal  frequencies,  at  which 
the  medium  gets  loaded,  are  given  by  the  modes’  absolute  values:  Ximpq  —  I  0Jimpti  I- 

The  positively-defined  forcing  frequency  Xjmpq  is  the  actual  physical  frequency  at  which 
the  Impq  term  in  the  expansion  for  the  tidal  potential  (or  stress  or  strain)  oscillates. 

The  motivation  for  keeping  also  the  modes  u>impq  is  subtle:  it  depends  upon  the  sign  of 
a>lmpq  whether  the  Impq  component  of  the  tide  lags  or  advances.  Specifically,  the  phase  lag 
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between  the  Impq  component  of  the  perturbed  primary’s  potential  U  and  the  hnpq  compo¬ 
nent  of  the  tide-raising  potential  W  generated  by  the  secondary  is  given  by 

^lmpq  ^ Impq  ^ Impq  I  ^ Impq  I  ^Hmpq  ^ lmpq  %-lmpq  ^ Impq  (,)hnpq  ’  (3^) 

where  the  time  lag  A  tlm  is  always  positive. 

While  the  lag  between  the  applied  stress  and  resulting  strain  in  a  sample  of  a  medium  is 
always  positive,  the  tidal  lag  can  be  positive  or  negative.  This  in  no  way  implies  violation  of 
causality  (the  time  lag  A tlm  is  always  positive).  Rather,  it  reflects  the  directional  difference 
between  the  planetocentric  positions  of  the  tide-raising  body  and  the  resulting  bulge.  For 
example,  the  principal  tide,  Impq  =  2200,  stays  behind  (has  a  positive  phase  lag  e-^oo)  when 
the  secondary  is  below  the  synchronous  orbit,  and  advances  (has  a  negative  phase  lag  e^oo ) 
when  the  secondary  is  at  a  higher  orbit.  To  summarise,  decomposition  of  a  tide  over  both 
positive  and  negative  modes  a>Im  (and  not  just  over  the  positive  frequencies  Ximpq )  does 
have  a  physical  meaning,  as  the  sign  of  a  mode  colmpq  carries  physical  information. 

Thus  we  arrive  at  the  following  conclusions: 

1.  As  the  fields  emerging  in  the  tidal  theory — the  tidal  potential,  stress,  and  strain — are  all 
real,  their  expansions  in  the  frequency  domain  may,  in  principle,  be  written  down  using 
the  positive  frequencies  /  only. 

2.  The  tidal  potential  contains  components  corresponding  to  the  tidal  modes  coIm  of  both 

the  positive  and  negative  signs.  While  the  Impq  components  of  the  potential,  stress,  and 
strain  oscillate  at  the  positive  frequencies  /;  =  \oJlmpq  |,  the  sign  of  a> y  points 

whether  the  lagging  of  the  Impq  component  of  the  bulge  is  positive  or  negative  (falling 
behind  or  advancing).  This  sign  enters  the  expression  for  the  appropriate  component  of 
the  torque  or  force.  Hence  a  consistent  tidal  theory  should  employ  expansions  over  both 
positive  and  negative  tidal  modes  u>lm  and  not  just  over  the  positive  Ximpq  ■ 

3.  In  order  to  rewrite  the  tidal  theory  in  terms  of  the  positively-defined  frequencies  Ximpq 
only,  one  must  inserts  “by  hand”  the  extra  multipliers 

sgn  oj,mpq  =  sgn  [(/  -2 p  +  q)n-  m  d]  (31) 

into  the  expressions  for  the  Impq  components  of  the  tidal  torque  and  force. 

4.  One  can  employ  a  rheological  law  (constitutive  equation  interconnecting  the  strain  and 
stress)  and  a  Navier-Stokes  equation  (the  second  law  of  Newton  for  an  element  of  the 
medium)  to  calculate  the  phase  lag  €Im  of  the  primary’s  potential  Ulm  relative  to  the 
potential  Wtm  generated  by  the  secondary.  If  both  these  equations  are  expanded,  in  the 
frequency  domain,  via  positively-defined  forcing  frequencies  Ximpq  only>  the  resulting 
phase  lag,  too,  will  emerge  as  a  function  of  Ximpq  '■ 

£ Impq  ^7  (%-lmpq  )*  (32) 

Within  this  treatment,  one  has  to  equip  the  lag,  “by  hand”,  with  the  sign  factor  (31).  On 
crossing  of  an  Impq  resonance,  the  factor  will  change  its  sign.  Accordingly,  the  Impq 
term  of  the  tidal  torque  (or  force)  will  change  its  sign  too. 

The  lag  (32)  is  the  negative  argument  of  the  complex  Love  number  kt  ( Ximpq  )■  Solution  of  the 
constitutive  and  Navier-Stokes  equations  renders  the  complex  Love  numbers,  from  which  one 
can  calculate  the  lags.  Hence  the  above  item  [4]  may  be  rephrased  in  the  following  manner: 

4'.  Under  the  convention  that  Ulmpq  =  U(x,mpq )  and  Wlmpq  =  W(xlmpq),  we  have: 

U Impq  =  kMimpq)  W tmpq  when  mlmpq  >  °’  i-e-when  colmpq  =  X,mpq .  (33a) 
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UlmPq  =  k*(Xlmpq)  Wlmpq  when  ulmpq<0,  i.e.,  when  colmpq  =  - Xlmpq ,  (33b) 

the  asterisk  denoting  the  complex  conjugation. 

This  ugly  convention,  a  switch  from  to  k*,  is  the  price  we  pay  for  employing  only  the 
positive  frequencies  in  our  expansions,  when  solving  the  constitutive  and  Navier-Stokes 
equations,  to  find  the  Love  number.  In  other  words,  this  is  a  price  for  our  pretending  that 
Wlm  and  Ulmpq  are  functions  of  Ximpq — whereas  in  reality  they  are  functions  of  u>lmpq . 

Alternative  to  this  would  be  expanding  the  stress,  strain,  and  the  potentials  over  the  positive 
and  negative  modes  colmpq ,  with  the  negative  frequencies  showing  up  in  the  equations.  With 
the  convention  that  Ulmpq  =  U(a>lmpq)  and  Wlmpq  =  W(wlmpq),  we  would  have 

Ulmpq  =kMmPq)Wlmpr  f()r  a11  W/mp9  •  (34) 


3.4  The  complex  rigidity  and  compliance:  stress-strain  relaxation 

The  stress  cannot  be  obtained  by  means  of  an  integral  operator  that  would  map  the  past 
history  of  the  strain,  U(f ')  over  t '  e  (— oo,  f],  to  the  value  of  §  at  time  t .  The  insufficiency 
of  such  an  operator  stems  from  the  presence  of  a  time-derivative  on  the  right-hand  side  of 
the  viscous  part  of  the  stress,  given  by  the  first  expression  in  (16).  Exceptional  are  the  cases 
of  no  viscosity  (e.g.,  a  purely  elastic  material). 

On  the  other  hand,  we  expect  on  physical  grounds  that  the  operator  J  inverse  to  /t  is  an 
integral  operator.  In  other  words,  we  assume  that  the  current  value  of  strain  depends  only  on 
the  present  and  past  values  taken  by  the  stress  and  not  on  the  current  rate  of  change  of  the 
stress.  The  assumption  works  for  weak  deformations,  i.e.,  insofar  as  no  plasticity  shows  up. 
So  we  assume  that  the  operator  J  is  just  an  integral  operator. 

Since  the  forced  medium  “remembers”  the  history  of  loading,  the  strain  at  time  t  must  be 
a  sum  of  small  installments  \  J (t  —  t')  dcryv(t each  of  which  stems  from  a  small  change 
d<jyV{t  —  r)  of  the  stress  at  an  earlier  time  t'  <  t.  The  entire  history  of  the  past  loading 
results,  at  the  time  t,  in  a  total  strain  uyv(t)  rendered  by  an  integral  operator  J {t)  acting  on 
the  entire  function  oyv(t ')  and  not  on  its  particular  value  (Karato  2008): 

OO  t 

2uyv(t)  =  J{t)  (Tyv  =  j  J( r)  <Tyv  (t  —  t)  dr  =  j  J(t  —  t')  oyv  (t1)  dt1,  (35) 

0  —  oo 

where  t 1  is  some  earlier  time  ( t '  <  t),  overdot  denotes  d/dt while  the  “age  variable” 
r  =  t  —  t  ’  is  reckoned  from  the  current  moment  t  and  is  aimed  back  into  the  past.  The 
so-defined  integral  operator  J(t)  is  called  the  compliance  operator,  while  its  kernel  J(t  —  t') 
goes  under  the  name  of  the  compliance  function  or  the  creep-response  function. 

Integrating  (35)  by  parts,  we  recast  the  compliance  operator  into  the  form  of 

oo 

2uyv(t)  =  J(t)oyv  =  J(0)(TyV{t)  -  J(oo)oyv(-oo)  +  ^  j  (t)  ffy v (t  -  r)dt  (36a) 

0 

t 

=  7(0)  ayv(t)  -  /(oo)  oyv(-oo)  +  j  j  (t  -  t')<TyV(t')dt'.  (36b) 

— OO 
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The  quantity  J  (oo)  is  the  relaxed  compliance.  Being  the  asymptotic  value  of  J  (t  —  t')  at 
t  —  t '  — >  oo,  this  parameter  corresponds  to  the  strain  after  complete  relaxation.  The  load  in 
the  infinite  past  may  be  assumed  zero,  and  the  term  —  J  (oo)  ayv(—o o)  may  be  dropped. 

The  second  important  quantity  emerging  in  (36)  is  the  unrelaxed  compliance  7(0),  which 
is  the  value  of  the  compliance  function  J(t  —  t')att  —  t'  =  0.  This  parameter  describes  the 
instantaneous  reaction  to  stressing,  and  thus  defines  the  elastic  part  of  the  deformation  (the 
rest  of  the  deformation  being  viscous  and  hereditary).  Thus  the  term  containing  the  unrelaxed 
compliance  7(0)  should  be  kept.  The  term,  though,  can  be  absorbed  into  the  integral  if  we 
agree  that  the  elastic  contribution  enters  the  compliance  function  not  as4 

J {t  —  t  ’)  =  7(0)  +  viscous  and  hereditary  terms,  (37) 

but  as 

J(t  —  t')  =  7(0)  ®(t  —  t ')  +  viscous  and  hereditary  terms,  (38) 

the  Heaviside  step-function  0(f  —  t')  being  unity  for  t  —  t  ’  >  0  and  zero  for  t  —  t ’  <  0.  As 
the  derivative  of  the  step-function  is  the  delta-function  S(t  —  t we  can  write  (36b)  as 

t 

2uyV{t)  =  j  7  (t  —  t  ’)  ctyV(t  ')dt ’ ,  with  J(t  —  t')  containing  7(0)  ®(t  —  t ').  (39) 

—  OO 

Formulae  (35-36)  and  (39)  are  but  different  expressions  for  the  compliance  operator  J  acting 
as 

2  Uyv  ~~  J  &yvm  (40) 

Inverse  to  the  compliance  operator  is  the  rigidity  operator  ft  defined  through 

Cfi/y—  2/it/yy.  (41) 

Generally,  jl  is  not  an  integral  operator,  but  is  an  integro-differential  operator.  So  it  cannot 
take  the  form  of  cryv{t)  =  2  J fi(t  —  t ')  uyv{t ')  dt However  it  can  be  written  as 

t 

<7yv(t)  =  2  J  p{t  -t')UyV(t')dt' ,  (42) 

—  OO 

if  we  permit  the  kernel  p(t  —  t')  to  contain  a  term  r/S(t  —  t'),  where  S(t  —  t')  is  the 
delta-function.  After  integration,  this  term  will  furnish  the  viscous  part  of  the  stress,  2  >7  uyv. 

The  kernel  p(t  —  t ')  is  termed  the  stress-relaxation  function.  Its  time-independent  part 
is  /x(0)  ®(t  —  t where  the  unrelaxed  rigidity  p(0)  is  inverse  to  the  unrelaxed  compliance 
7(0)  and  describes  the  elastic  part  of  deformation.  Each  term  in  p(t  —  t '),  which  neither  is 
a  constant  nor  contains  a  delta-function,  is  responsible  for  hereditary  reaction. 


4  Expressing  the  stress  through  the  strain,  we  encountered  three  possibilities:  the  elastic  stress  was  simply 
proportional  to  the  strain,  the  viscous  stress  was  proportional  to  the  time-derivative  of  the  strain,  while  the 
hereditary  stress  was  expressed  by  an  integral  operator  ft.  However,  when  we  express  the  strain  through  the 
stress,  we  place  the  viscosity  into  the  integral  operator,  so  the  viscous  reaction  also  looks  hereditary.  It  is  our 
convention,  though,  to  apply  the  term  hereditary  to  delayed  reactions  other  than  purely  viscous. 
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3.5  Stress-strain  relaxation  in  the  frequency  domain 

The  complex  compliance  J  (x)  and  the  complex  rigidity  Jx{x)  are,  by  definition,  the  Fourier 
images  not  of  the  J( r)  and  /r(r)  functions,  but  of  their  time-derivatives:5 


J(x)e'XTdx  =J  CO.  where  J (x)  =  /  i  (O  e  !J<r  dx. 


(43) 


and 


ti(x)e‘XTdx  =  A(r),  where  £i(x)  =  /  £(r)e  !;(rr/r, 


(44) 


o  o 

the  integrations  over  r  spanning  the  interval  [0,  oo),  as  both  kernels  are  nil  for  r  <  0  anyway. 
In  (43)  and  (44),  we  made  use  of  the  fact  (explained  in  Sect.  3.3)  that,  when  expanding  real 
fields,  it  is  sufficient  to  use  only  positive  frequencies. 

Formula  (35),  in  combination  with  the  Fourier  expansions  (23)  and  with  (43),  furnishes: 


'Xt 


>xt 


uyl,(x)e  dx  =  /  (Tuv(x)  J(X)e  dx 


(45) 


which  leads  us  to: 

2  uyv(x )  =  J(x)  dyy(x). 

Similarly,  insertion  of  (23)  into  (42)  leads  to  the  relation 

oyv(x)  =  2  /i(x)  Myy(x), 

comparison  whereof  with  (46)  immediately  entails: 

/(x)  A(x)  =  i- 

Writing  down  the  complex  rigidity  and  compliance  as 

£(x)  =  Im(x)I  exp|i'<5(x)| 

and 

J(X)  =  |/(X)I  exp[-i5(x)] , 
we  split  (48)  into  two  expressions: 


(46) 

(47) 

(48) 

(49) 

(50) 


]/(x)l 


1 

IA(x)l 


(51) 


and 


<pu(x)  =  Mx)  -  ^(x)- 


(52) 


From  the  latter,  we  see  that  the  angle  S  (x )  is  a  measure  of  lagging  of  a  strain  harmonic  mode 
relative  to  the  appropriate  harmonic  mode  of  the  stress.  It  ensues  from  (49-50)  that 


tan5(x) 


1m  [i(x)] 

Ke[J(x)] 


T^-e  [/x(x>] 


(53) 


5  Recall  that  it  is  the  time-derivative  of  J( r)  that  is  the  kernel  of  the  integral  operator  (39).  Hence,  to  arrive 


at  (46),  we  have  to  define  J(x)  as  the  Fourier  image  of  J  (r). 
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4  Complex  Love  numbers 


The  developments  presented  in  this  section  will  rest  on  a  theorem,  which  is  known  as  the 
correspondence  principle  or  the  elastic-viscoelastic  analogy,  and  which  applies  to  linear 
deformations  in  the  absence  of  nonconservative  (inertial)  forces.  While  the  literature  attri¬ 
butes  the  authorship  of  the  theorem  to  different  scholars,  its  true  pioneer  was  Sir  George 
Darwin  (1879).  One  of  the  corollaries  ensuing  from  this  theorem  is  that,  in  the  frequency 
domain,  the  complex  Love  numbers  are  expressed  via  the  complex  rigidity  or  compliance  in 
the  same  way  as  the  static  Love  numbers  are  expressed  via  the  relaxed  rigidity  or  compliance. 

As  was  pointed  out  by  Biot  (1954,  1958),  the  theorem  is  inapplicable  to  non-potential 
forces.  Hence  the  said  corollary  fails  in  the  case  of  librating  bodies,  because  of  the  presence 
of  the  inertial  force6  -aixrp,  where  p  is  the  density  and  « is  the  libration  angular  velocity.  So 
the  expression  (3)  for  the  Love  numbers,  generally,  cannot  be  employed  for  librating  bodies. 

Section  4. 1  below  explains  the  transition  from  the  stationary  Love  numbers  to  their  dynam¬ 
ical  counterparts,  the  Love  operators.  We  present  this  formalism  in  the  frequency  domain, 
in  the  spirit  of  Zahn  (1966)  who  pioneered  this  approach  in  application  to  a  purely  vis¬ 
cous  medium.  Section  4.2  addresses  the  negative  tidal  modes  emerging  in  the  Darwin-Kaula 
expansion  for  tides.  Employing  the  correspondence  principle,  in  Sect.  4.3  we  write  down  the 
expressions  for  the  factors  |k;(x)|  sine^x)  =  — Xnt[fc;(x)]  emerging  in  the  expansion  for 
tides.  Technical  details  of  this  derivation  are  discussed  in  Sects.  4.4-4. 5. 

For  more  on  the  correspondence  principle  and  its  applicability  to  Phobos  see  Appendix  B. 

4. 1  From  the  Love  numbers  to  the  Love  operators 


A  homogeneous  incompressible  primary,  when  perturbed  by  a  static  secondary,  yields  its 
form  and,  consequently,  has  its  potential  changed.  The  Ith  spherical  harmonic  Uj  (r)  of  the 
resulting  increment  of  the  primary’s  exterior  potential  is  related  to  the  Ith  spherical  harmonic 
W/(R,  r)  of  the  perturbing  exterior  potential  through  (2). 

Under  evolving  disturbances,  the  Love  numbers  become  operators  mapping  the  function 
W/{ R,  r*,  t '),  with  t '  e  (— oo,  r],  onto  a  value  of  U\  at  time  t: 

U,( r,  t )  =  1- j  k,(t)  Wj(R,  r*,  t’).  (54) 

Being  linear  for  weak  forcing,  the  operators  must  read: 


U,( r,  t)  = 


1+ 1  t_7.°° 

/  k/( r)  Wi  (R,  r*,  t  —  r)  dr 

r=0 

/+ 1  r'f 

/  kilt  —  t')  Wi  (R,  r* ,  t')dt' 

t  '=— oo 


(55a) 


or,  after  integration  by  parts: 


U,ir,  t ) 


\ki(())W(t) 


(R\l+l 

k/(oo)W(— oo)]+  I  —  J 


OO 

J  kjl r)  W]( R,  r*,  t  —  r)dr 
o 


(55b) 


6  The  centripetal  term  is  potential  and  causes  no  troubles,  except  for  the  need  for  a  degree-0  Love  number. 
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=  l-j  [kt(0)W(t) -£/(oo)W(-oo)]+(-j  /  k,(t  -  t')Wi(R,r*,t')dt' 

— OO 

(55c) 

/R\l+l  ( R\l+l  f  d  , 

=  “(7)  ki(oo)W(-oo)  +  j  -[kiff-t) 

—  OO 

-fc/(0)+  *i(O)0(f-f')]  W/(R,  r*,7) dt'.  (55d) 

Just  as  in  the  case  of  the  compliance  operator  (35-36),  in  expressions  (55)  we  obtain  the 
terms  ki(0)W(t)  and  —  kiipo)  W(— 00).  Of  the  latter  term,  we  can  get  rid  by  setting  W  (— 00) 
nil,  while  the  former  term  may  be  incorporated  into  the  kernel  in  exactly  the  same  way  as 
in  (37-39).  Thus,  dropping  the  unphysical  term  with  W (—00),  and  inserting  the  elastic  term 
into  the  Love  number  not  as  k/( 0)  but  as  k/( 0)  &(t  —  t we  simplify  (55d)  to 

Ui (r,  t)  =  (  7  )  k,(t -t')Wi(R,  r*,  t')dt',  (56) 

— OO 

with  k/(t  —  t ')  now  including,  as  its  part,  k/( 0)  &  (t  —  t')  instead  of  k/( 0). 

Were  the  body  perfectly  elastic,  k]{t  —  t ')  would  consist  of  the  instantaneous-reaction  term 
kj  (0)  &{t  —  t')  only.  Accordingly,  the  time-derivative  of  k\  would  be:  kj  (t  —  t')  =  kj  S(t  —  t ') 
where  k/  =  k/( 0),  so  expressions  (55-56)  would  coincide  with  (2). 

Similarly  to  introducing  the  complex  compliance,  one  can  define  the  complex  Love  num¬ 
bers  as  Fourier  transforms  of  kj  (t): 

CXD 

J  k^yxTdx  =ici  (r),  (57) 

0 


the  overdot  standing  for  d/dx .  Churkin  (1998)  suggested  to  term  the  time-derivatives  kj  ( t ) 
as  the  Love  functions  J  Inversion  of  (57)  trivially  yields: 

OO  OO 

k[(x)  =  j  ic,( r)  e~,XT  dz  =  fc;(0)  +  i  /  j  [^(r)  —  fy(0)  0(r)]  e~'xr  dr,  (58) 
0  0 

where  we  integrated  only  from  0  because  the  future  disturbance  contributes  nothing  to  the 

present  distortion,  so  k]{ r)  vanishes  at  r  <  0.  Recall  that  the  time  r  denotes  the  difference 
t  —  t' .  So  r  is  reckoned  from  the  present  moment  t  and  is  directed  back  into  the  past. 

Defining  in  the  standard  manner  the  Fourier  components  Ul  (x)  and  Wt  (x)  of  functions 
Uj  ( t )  and  Wj  (t),  we  write  (55)  in  the  frequency  domain: 

/R\l+1- 

C//CX)  =  (  — )  k,(x)  Wj(x).  (59) 

where  we  denote  the  frequency  simply  by  x  instead  of  the  awkward  Ximpq  ■  To  employ  (59) 

in  the  tidal  theory,  one  has  to  know  the  frequency-dependencies  kj(x)- 


7  Churkin  (1998)  used  functions  kj(t)  which  were,  due  to  a  difference  in  notations,  the  same  as  our  kj  (r). 
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4.2  The  positive  forcing  frequencies  y  =  M  versus  the  positive  and  negative  tidal  modes  a> 

It  should  be  remembered  that,  by  relying  on  formula  (59),  we  place  ourselves  on  thin  ice, 
because  the  similarity  of  this  formula  to  (46)  and  (47)  is  deceptive. 

In  (46)  and  (47),  it  was  legitimate  to  limit  our  expansions  of  the  stress  and  the  strain 
to  positive  frequencies  y  only.  Had  we  carried  out  those  expansions  over  both  positive  and 
negative  frequencies  u>,  we  would  have  obtained,  instead  of  (46)  and  (47),  similar  expressions 

2  uyv(a>)  =  J(a>)  dyv(a>)  and  ayv(co)  =  2  f2(a>)  Myv(a>).  (60) 

For  positive  a>,  these  would  simply  coincide  with  (46)  and  (47),  had  we  renamed  co  as  y .  For 
negative  a>  =  —y,  the  resulting  expressions  would  read  as 

2  Myv(-X)  =  H~X)  Vyv(-X)  and  dyv(-y)  =  2  A(~X)  uyv(-y),  (61) 

where  we  stick  to  the  agreement  that  y  always  stands  for  a  positive  quantity.  In  accordance 
with  (26),  complex  conjugation  of  (61)  would  then  return  us  to  (60). 

Physically,  the  negative-frequency  components  of  the  stress  or  strain  are  nonexistent.  If 
brought  into  consideration,  they  are  obliged  to  obey  (26)  and,  thus,  should  play  no  role, 
except  for  a  harmless  renormalisation  of  the  Fourier  components  in  (28). 

When  we  say  that  the  physically  measurable  stress  ayv  ( t )  is  equal  to  7 Ze  ^&yv  (x)  e  *  j , 

('/r 

it  is  unimportant  to  us  whether  the  y  -contribution  in  ayv  (r)  comes  from  the  term  cryv(y)  e 

only,  or  also  from  the  term  ayv(—y)  e  .  Indeed,  the  real  part  of  the  latter  is  a  clone  of 
the  real  part  of  the  former  (and  it  is  only  the  former  term  that  is  physical).  However,  things 
remain  that  simple  only  for  the  stress  and  the  strain. 

As  we  emphasised  in  Sect.  3.3,  the  situation  with  the  potentials  is  drastically  different. 

r  -  i'x?  I 

While  the  potential  U ( t )  is  still  equal  to  2^  T^-e  \U{y)e  ,  it  is  now  important  to  distin- 

iXt 

guish  whether  the  / -contribution  into  U (t)  comes  from  the  term  Uyv(y)  e  or  from  the 

i'(— x)r 

term  U (—y)  e  ,  or  from  both.  Although  the  negative  mode  (—  y)  would  bring  the  same 

input  as  the  positive  mode  /,  the  inputs  will  contribute  differently  into  the  tidal  torque. 
As  can  be  seen  from  (109),  the  secular  part  of  the  torque  is  proportional  to  sine/,  where 
ez  =  a>i  A timpq ,  with  the  time  lag  A tlmpq  being  positively  defined — see  formula  (105). 
So  the  secular  part  of  the  torque  explicitly  contains  the  sign  of  the  tidal  mode  a>impq  ■ 

Thus,  as  explained  in  Sect.  3.3,  a  more  accurate  form  of  formula  (59)  should  be: 

Ul(a>)=kl(a>)Wl(a>),  (62) 

where  u>  can  be  of  any  sign. 

If  however  we  pretend  that  the  potentials  depend  on  the  physical  frequency  y  =  \u>\  only, 
i.e.,  if  we  always  write  U (a>)  as  U(y),  then  (59)  must  be  written  as: 

U,(x)  =  W  W,(y ),  when  y  =  M  for  w  >  0,  (63a) 

and 

U,(y)  =k*(y)  W,(y),  when  x  =  M  for  a>  <  0.  (63b) 

Unless  we  keep  this  detail  in  mind,  we  shall  get  a  wrong  sign  for  the  Impq  component  of  the 
torque  after  the  despinning  secondary  crosses  the  appropriate  commensurability.  (We  shall, 
of  course,  be  able  to  mend  this  by  simply  inserting  the  sign  sgn  ojlmp/j  by  hand.) 
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4.3  The  complex  Love  number  as  a  function  of  the  complex  compliance 


While  the  static  Love  numbers  depend  on  the  static  rigidity  p  via  (3),  it  is  not  readily  appar¬ 
ent  that  the  same  relation  link  ki(x)  t°  d(x)-  The  correspondence  principle  discussed  in 
Appendix  B  tells  us  that  in  many  situations  the  viscoelastic  operational  moduli  /7(x),  J (x) 
obey  the  same  algebraic  relations  as  the  elastic  parameters  p,  J .  This  is  why  in  these  situa¬ 
tions  the  Fourier  or  Laplace  transform  of  our  viscoelastic  equations  mimic  (144),  except  that 
the  functions  acquire  overbars:  a  =  2puyv,  etc.  So  the  solution  is  Ut  =  kj  Wp  with  k\ 
retaining  the  same  dependence  on  p,  R  and  p  as  in  (3),  except  that  now  p  have  an  overbar: 


h(x) 


3 

2(1-  1) 
3 

2(1-  1) 


1 


1  ,  (2  Z~  +  4 1  +  3)  /x(x) 

+  IgpR 


1  _  3 
1  +  A,  J/Jix)  ~  2  (/  -  1) 


3  1 

2  (/  —  1)  1  +  Aip(x)/p 

JiX) 

Hx)  +  A,  J 


(64) 


Here  the  coefficients  A;  are  defined  via  the  unrelaxed  quantities  p  =  p( 0)  =  1/7  =  1/7(0) 
in  the  same  way  as  the  static  A;  were  introduced  via  the  static  (relaxed)  p  =  p(o o)  =  l/7(oo) 
in  (3). 

The  moral  of  the  story  is  that  at  low  frequencies  each  k\  depends  on  p  or  7  in  the  same 
way  as  the  static  k /  depends  on  the  static  p  or  7.  This  happens  because  at  low  frequencies 
we  drop  the  acceleration  term  in  the  equation  of  motion  (147b),  so  (147b)  looks  like  (144a). 

Representing  a  complex  Love  number  as 


ki(x)  =  Tie  [*Kx)]  +  i  Tm  [L/(x)] 


Mx)\  e 


-ieiix) 


(65) 


we  can  write  for  the  phase  lag  e/(x): 


tane/(x) 


or,  equivalently: 


1m  [ kt(x )] 
Tie  [£,(x)j 


\k(x)\  sine/(x)  =  -Tm  [A/(x)]- 


(66) 


(67) 


The  products  |£/(x)l  sine/(x)  standing  on  the  left-hand  side  in  (67)  emerge  also  in  the  Fou¬ 
rier  series  for  the  tidal  potential.  So  it  is  these  products,  not  k//Q,  that  enter  the  expansions 
for  forces,  torques  and  the  damping  rate.  They  link  the  body’s  rheology  with  its  spin  history: 
from  7(x)  to  ki(X )  t0  \ki(x)\  sine  (x),  the  latter  being  employed  in  the  theory  of  tides. 
Through  simple  algebra,  expressions  (64)  entail: 


|/q(x)|  sine/(x) 


-lm[ki(x)] 


_J1 _ —A/  7  Tm  [7 (x)] _ 

2  (Z  “  1}  {Tie  [7 (x)]  +  A/7)2  +  (Tm  [7_(x)])2  ' 

(68) 


As  we  know  from  Sects.  3.3  and  4.2,  formulae  (66-68)  should  be  used  with  care.  Since 
in  reality  U  and  k /  are  functions  not  of  x  but  of  a>,  then  formulae  (68)  should  be  equipped 
with  multipliers  sgn  u>lm  ,  when  plugged  into  the  expression  for  the  Impq  component  of 
the  tidal  force  or  torque.  This  prescription  is  equivalent  to  (63). 
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4.4  Should  we  write  kimpq  and  eimpq,  or  would  kj  and  6;  be  enough? 

In  the  preceding  subsection,  the  static  relation  (2)  was  generalised  to  evolving  settings  as 

U,mpq( r,  t)  =  (  -J  ki(t)  Wlmpq( R,  r*,  t'),  (69) 

where  Impq  is  a  quadruple  of  integers  employed  to  number  a  Fourier  mode  in  the  Darwin- 
Kaula  expansion  (96)  of  the  tide,  while  Uimpq  (r,  t)  and  W/mpq(R,  r*,  t ')  are  the  harmonics 
containing  cos (ximpqt  -  eimpq )  and  cos (ximpqt')  correspondingly. 

One  might  be  tempted  to  generalise  (2)  even  further  to 

/  R\ 

UlmpqiXi  t)  —  (  I  klmpq  (t)  Wlmpq  (R,  F  ,  t  ), 

with  the  Love  operator  (and,  consequently,  its  kernel,  the  Love  function)  bearing  dependence 
upon  m,  p,  and  q.  Accordingly,  (59)  would  become 

Ulmpq(x)  —  klmpq  (X)  Wlmpq  (x)*  (70) 

Fortunately,  insofar  as  the  Correspondence  Principle  is  valid,  the  functional  form  of  ki,npq  (x ) 
depends  upon  /  only  and,  thus,  can  be  written  down  simply  as  ki(ximpq)-  We  know  this  from 
the  considerations  offered  after  equations  (144).  There  we  explained  that  ki  depends  on 
X  =  Xlmpq  via  7(x),  while  the  functional  form  of  k /  bears  no  dependence  on  m,  p,  q. 

The  phase  lag  is  often  denoted  as  f-impq,  a  time-honoured  tradition  set  by  Kaula  (1964). 
However,  as  the  lag  is  expressed  through  ki  via  (66),  all  said  above  about  ki  applies  to  the 
lag:  while  the  functional  form  of  the  dependency  eimpq{ x)  may  be  different  for  different  Is, 
it  is  invariant  under  the  other  three  integers,  so  the  notation  e?  iximpq)  is  more  adequate. 

It  should  be  mentioned,  though,  that  for  bodies  of  pronounced  non-sphericity  coupling 
between  the  spherical  harmonics  furnishes  the  Love  numbers  and  lags  whose  expressions 
through  the  frequency,  for  a  fixed  Z,  have  different  functional  forms  for  different  m,  p,  q.  In 
these  cases,  the  notations  kunpq  and  eimpq  become  necessary  (Smith  1974;  Wahr  1981a, b,c; 
Dehant  1987a,b).  For  a  slightly  non-spherical  body,  the  Love  numbers  differ  from  the  Love 
numbers  of  the  spherical  reference  body  by  a  term  of  the  order  of  the  flattening,  so  a  small 
non-sphericity  can  usually  be  neglected. 


4.5  Rigidity  versus  self-gravitation 


For  small  bodies  and  small  terrestrial  planets,  the  values  of  Al  vary  from  about  unity  to 
hundreds.  For  example,  A 2  is  about  2  for  the  Earth,  20  for  Mars,  80  for  the  Moon,  and  200 
for  Iapetus.  For  superearths,  the  values  will  be  much  smaller  than  unity,  though.  Insofar  as 


one  can  approximate  (64)  with 


Ai  — -  »  1, 

I  7(X)  I 


(71) 


ki(x) 


3  J  (X) 

2(7  -  1)  7  (X)  +  A,  J 


3  J  (X) 
2  Ai  J 


+  0(|  J /{A]  i)|2). 


(72) 


except  in  the  closest  vicinity  of  an  Impq  resonance,  where  Ximpq  approaches  nil  and  J 
diverges  for  some  rheologies — like,  for  example,  for  those  of  Maxwell  or  Andrade. 
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Whenever  the  approximate  formula  (72)  is  applicable,  we  can  rewrite  (66)  as 


tane(x) 


1m  \ki(/_  )] 

Ke  [iq  (/ )] 


1m  [  7  (x  )] 

Ke  [7(x)j 


tanS(x), 


(73) 


wherefrom  we  readily  deduce  that  the  phase  lag  e(x)  of  the  tidal  frequency  x  coincides  with 
the  phase  lag  of  the  complex  compliance: 


e(x)  «  «(X).  (74) 

if  x  is  not  too  small  (so  we  are  not  too  close  to  the  commensurability).  This  way,  insofar  as 
the  condition  (71)  is  fulfilled,  the  component  f//(x)  of  the  primary’s  potential  lags  behind 
the  component  W/(x)  of  the  perturbed  potential  by  the  same  phase  angle  as  the  strain  lags 
behind  the  stress  at  frequency  x  in  a  sample  of  the  material.  Dependent  on  the  rheology,  a 
vanishing  tidal  frequency  may  or  may  not  limit  the  applicability  of  (71)  and  thus  cause  a 
considerable  difference  between  e  and  S. 

In  other  words,  the  approximation  is  valid  insofar  as  changes  of  shape  are  determined 
solely  by  the  local  material  properties,  and  not  by  self-gravitation  of  the  object  as  a  whole. 
This  depends  on  the  rheological  model.  For  a  Voigt  or  SAS8  solid  in  the  limit  of  x  — »•  0,  we 
have/(x)  — >■  /,  so  the  zero-frequency  limit  ofki(x)  is  the  static  Love  number/:/  =  |£(0)|.ln 
this  case,  approximation  (72-74)  remains  applicable  all  the  way  down  to  x  =0.  For  the  Max¬ 
well  and  Andrade  models,  however,  one  obtains,  for  vanishing  frequency:  7(x)  ~  l/0?x)> 
whence  /x  ~  r)x  and  ^(X )  approaches  the  hydrodynamical  Love  number  =  3/2. 

We  see  that  for  the  Voigt  and  SAS  models  approximation  (74)  works,  for  A/  q>>  1,  at 
all  frequencies,  as  the  condition  A/  1  can  be  set  for  any  value  of  x  •  For  the  Maxwell 
and  Andrade  solids,  the  condition  holds  only  at  x  >  rJ/1A; ^  A~  1 ,  and  so  does  the 
approximation  (74).  At  frequencies  below  this  threshold,  self-gravitation  “beats”  the  local 
material  properties,  and  the  behaviour  of  the  tidal  lag  deviates  from  that  of  the  lag  in  a  sample. 
This  fact  should  be  kept  in  mind  when  one  wants  to  explore  crossing  of  a  resonance. 

A  standard  caveat  is  in  order,  concerning  formulae  (72-74).  As  in  reality  the  potential  U 
is  a  function  of  a>  and  not  x ,  our  illegitimate  use  of  x  should  be  compensated  by  multiplying 
€l  (X-impq )  with  s8n  wimpq  >  when  e:  shows  up  in  the  expression  for  the  tidal  force  or  torque. 

4.6  The  case  of  inhomogeneous  bodies 

Tidal  dissipation  within  a  multilayer  near-spherical  body  is  studied  through  expanding  the 
involved  fields  over  the  spherical  harmonics  in  each  layer,  setting  the  boundary  conditions 
on  the  outer  surface,  and  using  the  matching  conditions  on  boundaries  between  layers.  This 
formalism  was  developed  by  Alterman  et  al.  (1959).  An  updated  discussion  of  the  method 
can  be  found  in  Sabadini  and  Vermeersen  (2004).  For  a  brief  review,  see  Legros  et  al.  (2006). 

Calculation  of  tidal  dissipation  in  a  Jovian  planet  is  an  even  more  formidable  task  (see 
Remus  et  al.  2012a  and  references  therein).  Flowever  dissipation  in  a  giant  planet  with  a  solid 
core  may  turn  out  to  be  approachable  by  analytic  means  (Remus  et  al.  2011,  2012b). 


The  acronym  SAS  stands  for  the  Standard  Anelastic  Solid,  also  known  as  the  Hohenemser-Prager  model. 
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5  Dissipation  at  different  frequencies 


5.1  The  data  collected  on  the  Earth:  in  the  lab,  over  seismological  basins,  and  through 
geodetic  measurements 


In  Efroimsky  and  Lainey  (2007),  we  considered  the  generic  rheological  model 

O' 

Q  =  (£x)  , 


(75a) 


where  x  is  the  tidal  frequency  and  £  is  a  parameter  having  the  dimensions  of  time.  The 
physical  meaning  of  this  parameter  is  elucidated  in  Efroimsky  and  Lainey  (2007).  Under  the 
special  choice  of  a  =  —  1  and  for  sufficiently  large  values  of  Q,  this  parameter  coincides  with 
the  time  lag  At  which,  for  this  special  rheology,  turns  out  to  be  the  same  at  all  frequencies. 

Actual  experiments  register  not  the  inverse  quality  factor  but  the  phase  lag  between  the 
reaction  and  the  action.  So  the  empirical  law  should  rather  be  written  down  as 


1 

sin  S 


Ol 

(£  x) 


(75b) 


which  is  equivalent  to  (75a),  provided  the  Q  factor  is  defined  there  as  Qenergy  and  not  as 
Q work — see  Sect.  2.2  for  details. 

The  applicability  realm  of  the  empirical  power  law  (75)  is  remarkably  broad — in  terms  of 
both  the  physical  constituency  of  the  bodies  and  their  chemical  composition.  Most  intriguing 
is  the  robust  universality  of  the  values  taken  by  the  index  a  for  very  different  materials: 
between  0.2  and  0.4  for  ices  and  silicates,  and  between  0.14  and  0.2  for  partial  melts. 
Historically,  two  communities  independently  converged  on  this  form  of  dependence. 

In  the  material  sciences,  the  rheological  model  (82),  wherefrom  the  power  law  (75b) 
stems,  traces  its  lineage  to  the  groundbreaking  work  by  Andrade  (1910)  who  explored  creep 
in  metals.  Through  the  subsequent  century,  this  law  was  found  to  be  applicable  to  a  vast 
variety  of  other  materials,  including  minerals  (Weertman  and  Weertman  1975;  Tan  et  al. 
1997)  and  their  partial  melts  (Fontaine  et  al.  2005).  With  almost  the  same  values  of  a,  this 
law  applies  also  to  ices  (McCarthy  et  al.  2007;  Castillo-Rogez  2009),  a  milestone  result, 
taken  the  physical  and  chemical  differences  between  ices  and  silicates.  It  is  agreed  upon 
that  in  crystalline  materials  the  Andrade  regime  can  find  its  microscopic  origin  both  in  the 
dynamics  of  dislocations  (Karato  and  Spetzler  1990)  and  in  the  grain-boundary  diffusional 
creep  (Gribb  and  Cooper  1998).  As  the  same  behaviour  is  inherent  in  metals,  silicates,  ices, 
and  even  glass-polyester  composites  (Nechada  et  al.  2005),  it  should  stem  from  a  single 
underlying  phenomenon  determined  by  some  principles  more  general  than  specific  material 
properties.  An  attempt  to  find  such  a  universal  mechanism  was  undertaken  by  Miguel  et  al. 
(2002).  See  also  the  theoretical  considerations  offered  in  Karato  and  Spetzler  (1990). 

In  seismology,  the  power  law  (75)  became  popular  in  the  second  part  of  the  XXth  century, 
with  the  progress  of  precise  measurements  on  large  seismological  basins  (Mitchell  1995; 
Stachnik  et  al.  2004;  Shito  et  al.  2004).  Further  confirmation  of  this  law  came  from  geodetic 
experiments  that  included:  (a)  satellite  laser  ranging  (SLR)  measurements  of  tidal  variations 
in  the  J2  component  of  the  gravity  field  of  the  Earth;  (b)  space-based  observations  of  tidal 
variations  in  the  Earth’s  rotation  rate;  and  (c)  space-based  measurements  of  the  Chandler 
wobble  period  and  damping  (Benjamin  et  al.  2006;  Eanes  and  Bettadpur  1996;  Eanes  1995). 
Not  surprisingly,  the  Andrade  law  became  a  key  element  in  the  recent  attempt  to  construct 
a  universal  rheological  model  of  the  Earth’s  mantle  (Birger  2007).  This  law  also  became  a 
component  of  the  non-hydrostatic-equilibrium  model  for  the  zonal  tides  in  an  inelastic  Earth 


Springer 


Bodily  tides  near  spin-orbit  resonances 


301 


by  Defraigne  and  Smits  (1999),  a  model  that  became  the  basis  for  the  IERS  Conventions 
(Petit  and  Luzum  2010).  While  the  lab  experiments  give  for  a  values  within  0. 2-0.4,  the 
geodetic  techniques  favour  the  interval  0.14-0.2.  This  minor  discrepancy  may  have  emerged 
due  to  the  presence  of  partial  melt  in  the  mantle  and,  possibly,  due  to  nonlinearity  at  high 
bounding  pressures  in  the  lower  mantle.  The  universality  of  the  Andrade  law  compels  us  to 
assume  that  (75)  is  generic  for  silicate  planets  and  moons.  Similarly,  the  applicability  of  (75) 
to  samples  of  ices  in  the  lab  is  likely  to  indicate  that  this  law  can  be  employed  for  description 
of  an  icy  moon  as  a  whole. 

Karato  and  Spetzler  (1990)  argue  that  at  frequencies  below  a  certain  threshold  xo  anelas- 
ticity  gives  way  to  purely  viscoelastic  behaviour,  so  the  parameter  a  becomes  close  to  unity. 
This  detail  is  missing  in  the  theory  of  the  Chandler  wobble  of  Mars,  by  Zharkov  and  Gudkova 
(2009).9  For  the  Earth’s  mantle,  the  threshold  corresponds  to  the  time-scale  about  a  year  or 
slightly  longer.  Although  in  Karato  and  Spetzler  (1990)  the  rheological  law  is  written  in  terms 
of  1/2,  we  shall  substitute  it  with  a  law  more  appropriate  to  the  studies  of  tides: 

-p 

kj  sine/=(£’x)  ,  where  p  =  0.2  —  0.4  for  x  >  Xo 

and  p  ~  1  for  x  <  Xo,  (76) 

X  being  the  frequency,  and  xo  being  the  frequency  threshold  below  which  viscosity  takes 
over  anelasticity. 

The  reason  why  we  write  the  power  scaling  law  as  (76)  and  not  as  (75)  is  that  at  the  lowest 
frequencies  the  geodetic  measurements  render  k\  sin  e/,  not  the  lag  angle  5  in  a  sample  (e.g., 
Benjamin  et  al.  2006).  For  the  same  reason,  we  denoted  the  exponents  in  (75)  and  (76)  with 
different  letters,  a  and  p.  Below  we  shall  see  that  these  exponents  do  not  always  coincide. 
Another  reason  for  giving  preference  to  (76)  is  that  not  only  the  sine  of  the  lag  but  also  the 
absolute  value  of  the  Love  number  is  frequency-dependent. 

5.2  Tidal  damping  in  the  Moon,  from  laser  ranging 

Fitting  of  the  LLR  data  to  the  scaling  law  (75)  was  carried  out  by  Williams  et  al.  (2001) 
who  demonstrated  that  the  lunar  mantle  possesses  quite  an  abnormal  value  of  the  exponent: 
—0.19.  A  later  reexamination  in  Williams  et  al.  (2008)  gave  a  less  embarrassing  value,  —0.09, 
which  was  still  negative  and  thus  seemed  to  contradict  our  knowledge  about  microphysical 
damping  mechanisms  in  minerals.  Thereupon,  Williams  and  Boggs  (2009)  commented: 

“There  is  a  weak  dependence  of  tidal  specific  dissipation  Q  on  period.  The  Q  increases 
from  ~30  at  a  month  to  ~35  at  one  year.  Q  for  rock  is  expected  to  have  a  weak  depen¬ 
dence  on  tidal  period,  but  it  is  expected  to  decrease  with  period  rather  than  increase. 
The  frequency  dependence  of  Q  deserves  further  attention  and  should  be  improved.” 

While  there  always  remains  a  possibility  of  the  raw  data  being  insufficient  or  of  the  fitting 
procedure  being  imperfect,  the  fact  is  that  the  negative  exponent  obtained  in  Williams  and 
Boggs  (2009)  does  not  necessarily  contradict  the  law  (75).  Indeed,  the  exponent  obtained 
through  LLR  was  not  the  a  from  (75)  but  was  the  p  from  (76).  The  distinction  is  critical  due 
to  the  difference  in  frequency-dependence  of  the  seismic  and  tidal  dissipation.  It  turns  out 
that  the  near-viscous  value  p  ~  I  from  the  second  line  of  (76),  appropriate  for  low  frequen¬ 
cies,  does  not  retain  its  value  all  the  way  to  the  zero  frequency.  Specifically,  in  Sect.  5.4  we 

9  This  circumstance  was  ignored  also  by  Defraigne  and  Smits  (1999).  Accordingly,  if  the  claims  by  Karato 
and  Spetzler  (1990)  are  correct,  the  table  of  corrections  for  the  tidal  variations  in  the  Earth’s  rotation  in  the 
IERS  Conventions  is  likely  to  contain  increasing  errors  for  periods  of  about  a  year  and  longer. 
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shall  see  that  below  the  frequency  l/(rM  At)  the  exponential  p  begins  to  decrease  with  the 
decrease  of  the  frequency.  (Here  xM  =  r\/ p  is  the  Maxwell  time,  rj  and  p  being  the  mantle’s 
viscosity  and  rigidity.)  As  the  frequency  becomes  lower,  p  changes  its  sign  and  eventually 
becomes  —  1  in  a  close  vicinity  of  x  =  0.  This  behaviour  follows  from  calculations  based  on 
a  realistic  rheology  (formulae  88-90  below),  and  goes  along  with  the  evident  physical  fact 
that  the  average  tidal  torque  must  vanish  in  a  resonance.10  In  Sect.  5.7,  comparison  of  this 
behaviour  with  the  LLR  results  will  yield  us  an  estimate  for  the  mean  lunar  viscosity. 

5.3  The  Andrade  model  as  an  example  of  viscoelastic  behaviour 

The  complex  compliance  of  a  Maxwell  material  contains  a  term  J  =  J( 0)  responsible  for 
the  elastic  part  of  the  deformation  and  a  term  —i/xv  describing  the  viscosity.  Any  additional 
term  in  the  compliance  will  correspond  to  other  forms  of  hereditary  reaction.  The  available 
geophysical  data  strongly  favour  a  particular  extension  of  the  Maxwell  approach,  the  And¬ 
rade  model  (Cottrell  and  Aytekin  1947;  Duval  1976).  In  modem  notations,  the  model  looks 


J (t  -  t ')  =  [j  +  ft  (t  -  t ')a  +  (t  -  t ')]  (77) 

where  a  is  a  dimensionless  parameter,  /3  is  a  dimensional  parameter,  ?;  is  the  steady-state 
viscosity,  while  J  is  the  unrelaxed  compliance  which  is  inverse  to  the  unrelaxed  rigidity: 
J  =  J(0)  =  1  / /x(0)  =  1  /  p.  So  (77)  is  the  Maxwell  model  amended  with  a  hereditary  term. 

An  example  illustrating  the  model  is  rendered  by  deformation  under  constant  loading.  The 
anelastic  term  dominates  at  short  times,  the  strain  being  a  convex  function  of  t  (the  so-called 
primary  or  transient  creep).  As  time  goes  on  and  the  load  is  kept  constant,  the  viscous  term 
becomes  larger,  and  the  strain  becomes  almost  linear  in  t  (the  secondary  creep). 

For  all  minerals  (including  ices),  a  assumes  values  from  0.14  to  0.4  (more  often,  through 
0.3) — see  the  references  in  Sect.  5.1.  The  parameter  /S  may  be  rewritten  as 

P  =  Jx-a  =  pTxx-a,  (78) 

the  quantity  xA  having  dimensions  of  time.  This  timescale  associated  with  the  Andrade  creep 
may  be  termed  as  the  “Andrade  time”  or  the  “anelastic  time”.  We  see  from  (78)  that  a  short 
xA  makes  the  anelasticity  more  pronounced,  while  a  long  xA  makes  the  anelasticity  weak.12 

Within  some  frequency  bands,  the  Andrade  time  gets  very  close  to  the  Maxwell  time 
(Castillo-Rogez  et  al.  2011;  Castillo-Rogez  and  Choukroun  2010): 

xA^xM^fi^JxMa  =  J1-01  r,-*  =  »r“,  (79) 

where  the  viscoelastic  timescale  (the  relaxation  Maxwell  time)  is  given  by 

xM  =  -  =  i?  J-  (80) 

/X 


10  For  example,  the  principal  tidal  torque  r j  =  ^2200  actin§  on  a  secondary  must  vanish  on  crossing  the 

synchronous  orbit.  This  indeed  happens,  as  the  exponential  p  becomes  —1  in  the  close  vicinity  of  X2200  = 

11  As  we  integrate  over  t  —  t '  e  [  0,  00  ),  the  terms  P{t  —  t  r)a  and  r]~ 1  (t  —  t can  do  without  the  Heaviside 
step-function  @(t  —  t We  remind  though  that  the  first  term,  J,  does  need  this  multiplier,  so  that  insertion  of 
(77)  into  (39)  renders  the  desired  J  8  (t  —  t')  under  the  integral,  after  the  differentiation  in  (39)  is  performed. 

12  While  the  Andrade  creep  is  caused  by  “unpinning”  of  jammed  dislocations  (Karato  and  Spetzler  1990; 
Miguel  et  al.  2002),  it  is  not  apparently  clear  if  rA  can  be  identified  with  the  typical  time  of  unpinning. 
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One  cannot  expect  xA  and  xM  to  coincide  in  all  situations,  as  they  may  possess  some  degree  of 
frequency-dependence.  There  exist  indications  that  in  the  Earth’s  mantle  the  role  of  anelas- 
ticity  (compared  to  viscoelasticity)  undergoes  a  decrease  when  the  frequencies  become  lower 
than  1  /year — see  the  miscrophysical  model  suggested  in  section  5.2.3  of  Karato  and  Spetzler 
(1990).  The  relation  between  xA  and  xM  may  depend  also  upon  the  intensity  of  loading,  i.e., 
upon  the  damping  mechanisms  involved.  The  microphysical  model  considered  in  Karato  and 
Spetzler  (1990)  was  applicable  to  strong  deformations,  with  anelastic  dissipation  being  dom¬ 
inated  by  dislocations  unpinning.  Accordingly,  the  dominance  of  viscosity  over  anelasticity 
(xA  <3C  xM )  at  low  frequencies  may  be  regarded  proven  for  strong  deformations  only.  At  low 
stresses,  when  the  grain-boundary  diffusion  mechanism  is  dominant,  the  values  of  xA  and  xM 
may  remain  comparable  at  low  frequencies.  The  topic  needs  further  research. 

In  terms  of  the  Andrade  and  Maxwell  times,  the  compliance  becomes: 


J(t-t')=  J 


t  ')■ 


In  the  frequency  domain,  compliance  (81)  will  look: 


(81) 


~J(X)  =  J  +  PUxra  Td+or)-  — 

vx 

=  J  [1  +  (i  X  Tx)-“  r  (1  +  a)  -  i  (X  r Mr*]  , 


(82a) 

(82b) 


where  x  is  the  frequency  and  T  is  the  Gamma  function.  The  imaginary  and  real  parts  are: 

lm\J(x)]  =  ~(ri  x)~'  ~  X~a  P  sin  (a  jr /2)  T(a  +  1)  (83a) 

=  -J(X*mT1  ~  J  (X^a)-*  sin  ( a  7r /2 )  r(a  +  l),  (83b) 


TZe[J(X)]  =  J  +  X~“  P  cos  (a  tx/2)  T(a  +  1) 

=  J  +  J  (xTt)~“  cos  (a  jr/2)  T(a  +  1), 
whence  we  obtain  the  following  dependence  of  the  phase  lag  upon  the  frequency: 


tan<5(x)  =  - 


lm[j ix)\  _  tn  x)  +  X  yS  sin(a  ^/2)  r  (a  +  1) 


Tie  [7(x)] 


/z  +  x  P  cos  (“  rr/2)  T  (a  +  1) 


(84a) 

(84b) 


Z  1  C  +  z  “  sin  (a  tv /2)  T  (a  +  1) 

1  +  z”a  cos  (a  7 r /2)  T  (a  +  1 ) 

Here  z  is  the  dimensionless  frequency  defined  as 

z  =  X  rA  =  x  Tu  C, 

while  f  is  a  dimensionless  parameter  of  the  Andrade  model: 

f  =  xaIxm- 

5.4  Tidal  response  of  viscoelastic  spherical  bodies  obeying  the  Andrade  and  Maxwell 
models 


(85a) 

(85b) 

(86) 

(87) 


An  Impq  term  of  the  tidal  torque  is  proportional  to  k:(x)  sine;(x)  =  l^/(X/mp9)l  sin 
6/  (Ximpq )•  These  factors’  frequency  dependence  can  be  obtained  by  combining  (68)  with 
(82).  Leaving  details  for  Appendix  C.2,  we  present  the  results,  without  the  sign  multiplier. 
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•  The  high-frequency  band:  x  Xhi  > 

\k,(X)\  sine,(x)  ~  2(/3_  |}  +  1)2  sin  (°^)  T(a  +  \)  i;~a  (zM  x)““  •  (88) 

At  high  frequencies,  anelasticity  dominates.  So,  dependent  upon  the  microphysics  of  the 
mantle,  the  parameter  f  may  be  of  order  unity  or  slightly  lower.  We  say  slightly,  because  we 
expect  both  anelasticity  and  viscosity  to  be  present  near  the  transitional  zone.  (A  too  low  f 
would  eliminate  viscosity  from  the  picture.)  Were  the  parameter  f  evolving  slowly  with  the 
decrease  of  the  frequency,  we  would  be  able  to  say  that  for  small  bodies  and  small  terrestrial 
planets  (i.e.,  for  A;  1),  the  boundary  between  the  high  and  intermediate  frequencies  is 

a 

-1  „  l-« 

Xhi  —  rM  ?  > 

while  for  large  terrestrial  planets  (i.e.,  for  At  <£.  1)  the  boundary  is 

xH,  =  r~l  =  T~lrl- 

Both  these  formulae  are  derived  in  the  Appendix. 

In  reality,  these  two  estimates  may  be  overruled  by  microphysics.  At  a  certain  threshold 
frequency,  the  anelastic  timescale  rA  and  therefore  the  parameter  f  begin  to  grow  rapidly  with 
the  decrease  of  the  frequency.  This  happens  when  the  effectiveness  of  anelastic  dissipation 
(defect  unpinning)  begins  to  decrease.  Whenever  this  threshold  happens  to  be  higher  than  the 
above  estimates  for  xHi  ,  one  has  simply  to  identify  xHi  whh  the  said  threshold.  According 
to  Karato  and  Spetzler  (1990),  for  the  Earth’s  mantle  this  threshold  is  about  1/year.  Hence 
for  terrestrial  planets  it  may  be  more  realistic  to  set,  by  hand: 

Xhi  ~  1  /year- 

•  The  intermediate-frequency  band:  xHi  ^  X  ^  (A;  +  l)-1, 

IMX)I  sine;(X)  -  {Tm  x)-‘  .  (89) 

a 

While  Appendix  C.2  renders  z^1  t;  “  for  the  upper  bound,  here  we  identify  this  bound 
with  the  afore-mentioned  threshold  xHi  >  in  understanding  that  below  this  frequency  f  begins 
to  grow  rapidly,  thus  causing  a  switch  from  the  Andrade  to  a  near-Maxwell  regime. 

•  The  low-frequency  band:  z~l  (At  +  1)~*  x, 

\k,(x)\  sin e,(x)  «  2(/^_  r M  X-  (90) 

Scaling  laws  (88)  and  (89)  mimic,  up  to  constant  factors,  the  frequency-dependencies  of 
|/(X)|  sin<5(x)  =  —  Tm[J(x)]  at  high  and  low  frequencies,  correspondingly,  see  Eq.  (83). 

Expression  (90)  however  shows  a  remarkable  phenomenon  inherent  only  in  the  tidal  lag¬ 
ging,  and  not  in  the  lagging  in  a  sample  of  material.  The  |k/(x)l  sine;(x)  factor  changes 
its  behaviour  drastically  on  close  approach  to  the  zero  frequency.  Having  reached  a  finite 
maximum  at  about  x  =  t  ~ 1  (A[  +  1 )  1 ,  the  said  factor  begins  to  scale  linearly  in  x  as  x 
approaches  zero.  This  way,  \kt  (x )  I  sin  et  (x)  decreases  continuously  on  close  approach  to  a 
resonance  and  goes  through  nil  together  with  the  frequency.  So  neither  the  tidal  torque  nor 
the  tidal  force  explodes  in  resonances.  In  a  somewhat  heuristic  manner,  this  change  in  the 
frequency-dependence  was  pointed  out,  for  /  =  2,  in  Efroimsky  and  Williams  (2009). 
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Fig.  1  Tidal  response  of  a  homogeneous  spherical  Andrade  body,  set  against  the  decadic  logarithm  of  the 
forcing  frequency  y  (in  Hz).  The  blue  curve  renders  the  decadic  logarithm  of  the  absolute  value  of  the  quad- 
rupole  complex  Love  number,  ]g  ki  =  lg  Of/ 1 1  -  The  green  and  red  curves  depict  the  logarithms  of  the  real 
and  the  negative  imaginary  parts  of  the  Love  number:  lg  1Ze[k2ix )]  =  lg(&2  cos  f2)  and  lg{ — XmfeCx)])  = 
lg(— ^2  sin  €2),  accordingly.  The  change  in  the  slope  of  the  red  curve  (the  "elbow”),  which  takes  place  to  the 
right  of  the  maximum,  corresponds  to  the  switch  from  viscosity  dominance  at  lower  frequencies  to  anelasticity 
dominance  at  higher  frequencies.  The  parameters  A2  and  rM  were  given  values  appropriate  to  a  homogeneous 
Moon  with  a  low  viscosity,  as  described  in  Sect.  5.7.  The  plots  were  generated  for  an  Andrade  body  with 
0  =  1  at  all  frequencies.  Setting  the  body  Maxwell  at  lower  frequencies  will  only  slightly  change  the  shape 
of  the  “elbow”  and  will  have  virtually  no  effect  on  the  maximum 


In  that  paper,  we  also  tried  to  describe,  through  semi-qualitative  reasoning,  the  frequency- 
dependence  of  the  Love  number  ki.  We  concluded  that  it  should  be  proportional  to  the  cosine 
of  the  phase  lag.  Accurate  treatment  shows  that  for  the  Kelvin- Voigt  model  this  is  true  over 
all  frequencies.  For  the  Maxwell  model,  this  relation  works  in  the  low-frequency  limit.  For 
details,  see  the  Electronic  Supplementary  Material  5  “The  behaviour  of  kj(x)  =  |£/(x)|  in 
the  limit  of  vanishing  tidal  frequency  / ,  within  the  Andrade  and  Maxwell  models”. 


5.5  Example 

Figure  1  shows  the  absolute  value,  k2  =  |&2(x)l-  the  real  part,  ki  cos  e 2  =  Tie  [^2 (X )] ^  and 
the  negative  imaginary  part,  £2  sin  €2  =  —Tin  [^2 (x )] ^  °f  the  quadrupole  Love  number.  Each 
quantity’s  decadic  logarithm  is  plotted  against  the  decadic  logarithm  of  the  frequency  x  (in 
Hz).  The  curves  were  obtained  by  insertion  of  formulae  (83-84)  into  (64).  As  an  example, 
the  case  of  —Tin  [^2  (X )]  is  worked  out  in  Appendix  C.2,  formulae  (168-170). 

Both  in  the  high-  and  low-frequency  limits,  the  negative  imaginary  part  of  k2(x)  (the  red 
curve  in  Fig.  1)  approaches  zero.  Accordingly,  in  the  low-  and  high-frequency  bands  the  real 
part  (the  green  curve)  virtually  coincides  with  the  absolute  value  (the  blue  curve). 

While  on  the  left  and  on  the  close  right  of  the  peak,  dissipation  is  mainly  due  to  viscosity, 
friction  at  higher  frequencies  is  mainly  due  to  anelasticity.  This  switch  corresponds  to  the 
change  of  the  slope  of  the  red  curve  at  high  frequencies  (for  our  choice  of  parameters,  at 
around  1CV5  Hz).  This  change  of  the  slope  is  often  called  the  elbow. 

Figure  1  was  generated  for  A2  =  80.5andrM  =  3.75x  105  s.  The  value  of  A2  corresponds 
to  the  Moon  modeled  by  a  homogeneous  sphere  of  rigidity  /r  =  0.8  x  101 1  Pa.  Our  choice  of 
the  value  of  rM  =  rj/fz  corresponds  to  a  homogeneous  Moon  with  the  said  value  of  rigidity 
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and  with  viscosity  set  to  be  r]  =  3.75  x  105  Pa  s.  The  reason  why  we  consider  an  example 
with  such  a  low  value  of  q  will  be  explained  in  Sect.  5.7.  Finally,  it  was  assumed  for  simplicity 
that  f  =  1,  i.e.,  that  rA  =  xM.  Although  unphysical  at  low  frequencies,  this  simplification 
only  slightly  changes  the  shape  of  the  “elbow”  and  exerts  virtually  no  influence  upon  the 
maximum  of  the  red  curve,  provided  the  maximum  is  located  well  into  the  viscosity  zone. 

5.6  Crossing  a  resonance — with  a  chance  for  entrapment 

In  the  expansion  for  the  tidal  torque,  factors  (88-90)  appear  in  the  company  of  sgn  u>.  For 
example,  the  factor  (90)  enters  the  hnpq  term  of  the  torque  as 

I  (^linpq  )  I  sin  6/  (Xlmpq  )  ^Impq 

3  3 

2 |  )  ^  Tir  S£n  (*J[mpq  't  (J  —  1  )  ^  Tlf  '  (91) 

Naturally,  this  term  depends  on  colmpq,  not  on  xhnpq  =  \“>impq\-  Together  with  wlmpq,  the 
term  goes  continuously  through  zero,  and  changes  its  sign  as  the  linpq  resonance  is  crossed. 
This  gets  along  with  the  physically  evident  fact  that  the  principal  (i.e.,  2200)  term  of  the  tidal 
torque  should  vanish  as  the  secondary  approaches  the  synchronous  orbit. 

An  Impq  term  of  the  torque  changes  its  sign  and  thus  creates  a  chance  for  entrapment.  As 
the  value  of  an  Impq  term  of  the  torque  is  much  lower  than  that  of  the  principal,  2200  term, 
a  perfectly  spherical  body  will  never  get  stuck  in  a  resonance  other  than  2200.  (The  latter  is, 
of  course,  the  1  :  1  resonance  in  which  the  principal  term  of  the  torque  vanishes.)  However, 
the  presence  of  the  triaxiality-generated  torque  is  known  to  contribute  to  the  probabilities 
of  entrapment  into  other  resonances  (provided  the  eccentricity  is  not  zero).  Typically,  in  the 
literature  they  consider  a  superposition  of  the  triaxiality-generated  torque  with  the  principal 
tidal  term.  We  would  point  out  that  the  “trap”  shape  of  the  Impq  term  (91)  makes  this  term 
relevant  for  the  study  of  entrapment  in  the  Impq  resonance.  In  some  situations,  one  has  to 
take  into  account  also  the  non-principal  terms  of  the  tidal  torque. 

5.7  Comparison  with  the  LLR  results 

As  we  mentioned  above,  fitting  of  the  LLR  data  to  the  power  law  has  resulted  in  a  very 
small  negative  exponential  p  =  —0.19  (Williams  et  al.  2001).  As  the  measurements  of  the 
lunar  damping  described  in  Williams  et  al.  (2001)  rendered  information  on  the  tidal  and  not 
seismic  dissipation,  those  results  can  and  should  be  compared  to  the  scaling  law  (88-90).  As 
the  small  negative  exponential  was  devised  from  observations  over  periods  of  a  month  to  a 
year,  it  is  natural  to  presume  that  the  appropriate  frequencies  were  close  to  or  slightly  below 
the  frequency -  '  . .  at  which  the  factor  k 2  sin  has  its  peak,  as  in  Fig.  1: 

rM  !a2  + 

3  x  106  s  «  0.1  year  =  xM  (A2  +  1)  «  rM  A2  =  ^  A2  =  R)1  •  (92> 

Hence,  were  the  Moon  a  uniform  viscoelastic  body,  its  viscosity  would  be  only 

ij  =  3x  1016Pa  s.  (93) 

Hence  the  lower  lunar  mantle  contains  a  high  percentage  of  partial  melt,  a  fact  which  goes 
along  with  the  model  by  Weber  et  al.  (2011),  and  which  was  anticipated  in  Williams  et  al. 
(2001)  and  Williams  and  Boggs  (2009)  following  an  earlier  work  by  Nakamura  et  al.  (1974). 
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6  The  polar  component  of  the  tidal  torque  acting  on  the  primary 

Let  vector  r  =  (r,  k,  <f>)  point  from  the  centre  of  the  primary  toward  a  point-like  secondary 
of  mass  Msec.  Associating  the  coordinate  system  with  the  primary,  we  reckon  the  latitude 
4>  from  the  equator.  Setting  the  coordinate  system  to  corotate,  we  determine  the  longitude  k 
from  a  fixed  meridian.  The  tidally  induced  component  of  the  primary’s  potential,  U ,  can  be 
generated  either  by  this  secondary  itself  or  by  some  other  secondary  of  mass  M*ec  located  at 
r*  =  (r*,  A.*,  </>*).  In  either  situation,  the  tidally  induced  potential  U  generates  a  tidal  force 
and  a  tidal  torque  wherewith  the  secondary  of  mass  Msec  acts  on  the  primary. 

The  scope  of  this  paper  is  limited  to  low  i.  When  the  role  of  the  primary  is  played  by  a 
planet,  the  secondary  being  its  satellite,  i  is  the  satellite’s  inclination.  When  the  role  of  the 
primary  is  played  by  the  satellite,  the  planet  acting  as  its  secondary,  i  denotes  the  satellite’s 
obliquity.  Similarly,  when  the  planet  is  regarded  as  a  primary  and  its  host  star  is  treated  as  its 
secondary,  i  is  the  obliquity  of  the  planet.  In  all  these  cases,  the  smallness  of  i  indicates  that  the 
tidal  torque  acting  on  the  primary  is  reduced  to  its  polar  component.  The  other  components 
of  the  torque  will  be  neglected  in  this  approximation. 

The  polar  component  of  the  torque  acting  on  the  primary  is  the  negative  of  the  partial 
derivative  of  the  tidal  potential,  with  respect  to  the  primary’s  sidereal  angle: 

dU(r) 

T(r)  =  -Msec  y,  (94) 

ou 

9  standing  for  the  primary’s  sidereal  angle.  The  formula  is  convenient  when  the  tidal  potential 
U  is  expressed  through  the  secondary's  orbital  elements  and  the  primary’s  sidereal  angle. 1 

Here  and  hereafter  we  deliberately  refer  to  a  primary  and  a  secondary  in  lieu  of  a  planet 
and  a  satellite.  The  preference  stems  from  our  intention  to  extend  the  formalism  to  setting 
where  a  moon  is  playing  the  role  of  a  tidally-perturbed  primary,  the  planet  being  its  tide-pro¬ 
ducing  secondary.  Similarly,  when  addressing  the  rotation  of  Mercury,  we  interpret  the  Sun 
as  a  secondary  that  is  causing  a  tide  on  the  effectively  primary  body,  Mercury. 


7  The  tidal  potential 


7.1  Darwin  (1879)  and  Kaula  (1964) 


The  potential  produced  at  point  R=  (R,  k,  <j>)  by  a  secondary  body  of  mass  M*,  located  at 
r*  =  (r*,  k* ,  <p*)  with  r*  >  R,  is  given  by  (1).  When  a  tide-raising  secondary  located  at 
r*  distorts  the  shape  of  the  primary,  the  potential  generated  by  the  primary  at  some  exterior 
point  r  acquires  an  amendment  given  by  (2).  Insertion  of  (1)  into  (2)  entails 


13  Were  the  potential  written  down  in  the  spherical  coordinates  associated  with  the  primary’s  equator  and 
corotating  with  the  primary,  the  polar  component  of  the  tidal  torque  could  be  calculated  with  aid  of  the 
expression 


T  (r)  =  Msec 


dUjr) 
d  k 


derived,  for  example,  in  Williams  and  Efroimsky  (2012).  That  the  expression  agrees  with  (94)  can  be  seen 
from  the  formula 


k  =  —  9  +  £}  +  co  +  v+  0(i2)  =  —0  +  Q  +  ai  +  A4+2e  sinAf  4-  0(e2)  +  0(i2), 


e,  (,  <u,  £2,  v  and  .VI  being  the  eccentricity,  inclination,  argument  of  the  pericentre,  longitude  of  the  node,  true 
anomaly,  and  mean  anomaly  of  the  tide-raising  secondary. 
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vn.-GMuXffpaZ.fcsji 

1=2  i  '  m= 0 

X  (2  -  <$0m ) Pim (sin </>) Pim (sin  0*)  cos m(k  -  A*). 


21+1  I 

k/  R  (/  —  m)\ 


(95) 


A  more  practical  formula  was  offered  by  Kaula  (1961,  1964),  who  switched  from  the  spheri¬ 
cal  coordinates  to  the  Kepler  elements  (a*,  e* ,  i* ,  £2*,  co*,  ,M*)and(a,  e,  i,  £2,  co,  M) 
of  the  secondaries  located  at  r*  and  r.  Application  of  this  technique  to  (95)  results  in 


U (r)  =  - 


OO  J+l  J  1  1 

^  ( R\  GM*  /R\  +  ( l-m)\  + 

^  k  (a)  W  ^  (l  +  m)\  (2~S,jm)  XF^P{l  ) 

1=2  X  7  x  7  m= 0  v  7  p=0 


LXJ  l  LXJ 

X  Glpq(e*)X^mk(i)  X  cos[(^  -  m9*^j  -  ( vi,nhj  -  i h0)] 


q=—o o 


/z=0 


y=-oo 


(96) 


where 


^  =  (/  -  2/?)*;*  +  (/  -  2p  +  q)M*  +  m  £2*,  (97) 

u/m/y  =  (l  —  2 h)to  +  (/  —  2/7  +  y')A4  +  m  £2,  (98) 

9  =  6*  being  the  sidereal  angle,  Gipq  (e)  signifying  the  eccentricity  functions,14  and  Fimp(i) 
denoting  the  inclination  functions  (Gooding  and  Wagner  2008). 

Being  equivalent  for  a  planet  with  an  instant  response  of  the  shape,  (95)  and  (96)  disagree 
when  friction-caused  delays  come  into  play.  Kaula’s  expression  (96),  as  well  as  its  truncated, 
Darwin’s,  version,15  is  capable  of  accommodating  separate  phase  lags  for  each  mode: 


U (r)  =  - 


no  /+1  /  /  / 

/ R\  GM*  /  R\  +  (Z  —  m)\  + 

J  k  (  a  )  t )  Jo  ^+m)\  (2  “  S0m)  J  Flmp(l  } 


OO  l  o o 

^  "  Glpq(&  )  ^  '  Fimh  (l)  ^  '  Glhj(e)  COS  ^  {yimhj  WlO')  £/ra/?gji 


q= — 00 


h= 0 


j=-oo 


(99) 


where 

=  [(/  -  2p)  (b*  +  (1  -  2p  +  q)  M*  +  m  (£2*  -  0*)  ]  Ar/mp<? 

=  ^Impq  At Impq  =  ±  Xlmpq  Atlmpq  (100) 

is  the  phase  lag.  The  tidal  mode  (+*mpq  introduced  in  (100)  is 

tfmpq  =  (f-  2 p)  co*  +  (/  -  2p  +  * )  A1*  +  m  (6*  -  0*),  (101) 

while  the  positively-defined  quantity 


14  .  ...  .  (— /— 1),  (l—2p) 

Functions  Gihj(e )  coincide  with  the  Hansen  polynomials  ^  ^  (£).  In  the  Electronic  Supple¬ 

mentary  Material  2  “The  eccentricity  functions",  we  present  a  table  of  the  Gihj  {e)  functions,  which  are  needed 
for  expanding  the  tides  up  to  terms  with  e 6,  inclusively. 

While  the  treatment  by  Kaula  (1964)  entails  the  infinite  Fourier  series  (96),  the  development  by  Darwin 
(1879)  renders  its  partial  sum  with  |/|,  \q\,  \j\  <  2.  For  introduction  into  Darwin’s  method  see  Ferraz-Mello 
et  al.  (2008).  Be  mindful  that  these  authors’  convention  on  the  notations  r  and  r*  is  opposite  to  ours. 
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Ximpq  =  I  tfmpq  I  =  I  (Z  -  Ip)  <»*  +  (/  -  2p  +  9)  .M*  +  m  (£2*  -  0*)  |  (102) 

is  the  actual  physical  Impq  frequency  excited  by  the  tide  in  the  primary.  The  corresponding 
positively-defined  time  delay  At/mpq  =  Ati  (ximpq)  depends  on  this  physical  frequency,  the 
functional  forms  of  this  dependence  being  different  for  different  materials. 

In  neglect  of  the  apsidal  and  nodal  precession  and  also  of  Mo,  the  above  formulae  look: 


(0 Impq  =  (l  -  2p  +  q)  n  -  m  e, 

(103) 

Ximpq  =  1  COlmpq  \  =  \  (l  -  2p  +  q)  n  -  m  0  \, 

(104) 

^Impq 

—  t Olmpq  Atjmpq  ~  (/  2 p  A  q)  n  m  h  j  A t\mpq 

(105a) 

=  Ximpq  At,  (ximpq)  sgn  [(/  -2  p  +  q)n-md\. 

(105b) 

Formulae  (96)  and  (99)  constitute  the  pivotal  result  of  Kaula’s  theory  of  tides.  Importantly, 
this  theory  imposes  no  a  priori  constraint  on  the  form  of  frequency-dependence  of  the  lags. 

8  The  Darwin  torque 

The  popular  model  of  entrapment  into  the  1:1  resonance,  developed  by  Goldreich  (1966), 
rests  on  the  MacDonald  torque  and  thus  implies  an  unphysical  rheology  of  the  satellite’s 
material  (Williams  and  Efroimsky  2012).  The  rheology  is  given  by  (75)  with  a  =  —  1.  More 
realistic  is  the  dissipation  law  (76).  An  even  more  accurate  and  practical  formulation  of 
the  damping  law,  stemming  from  the  Andrade  formula  for  the  compliance,  is  rendered  by 
(88-90).  These  formulae  should  be  inserted  into  the  Darwin-Kaula  theory  of  tides. 


8. 1  The  secular  and  oscillating  parts  of  the  Darwin  torque 

8.1.1  The  general  formula 

Direct  differentiation  of  (99)  with  respect  to  —6  will  result  in  the  expression16 

'+'  .  .  l  i  .  i 


T  =  - 


Z(-)  (4)  i.  jrfZ  2"  Z  MO 

a*  V«V^oa  +  m)! 

oo  l  oo 

"  Glpq^C  )  y  'j  Ffrnh  (l)  y  'J  G7/jy(£)^/Sin^l^m^g  Vlmhj  ^Impq  J  •  (106) 


<7=— o o 


J=-o o 


If  the  tidally-perturbed  and  tide-raising  secondaries  are  the  same  body,  then  Msec  =  M*ec, 
and  all  the  elements  coincide  with  their  counterparts  with  an  asterisk.  Hence  the  differences 

Vlmpq  ~  vlmhj 

=  (/  —  2p  +  q)  M*  —  (l  —  2 h  +  j)  M  +  m  (£2*  —  £2)  +  l  (co*  —  at)  —  2  pa>*  +  2  hco 

(107) 

get  simplified  to 

v*mpq  -  vimhj  =  (2  h  -  2  p  +  q  -  j)  M*  +  (2  h  -  2  p)  co* ,  (108) 


16  For  justification  of  this  operation,  see  Sect.  6  in  Efroimsky  and  Williams  (2009). 
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an  expression  containing  both  short-period  contributions  proportional  to  the  mean  anomaly, 
and  long-period  contributions  proportional  to  the  argument  of  the  pericentre. 


8.1.2  The  secular,  the  purely-short-period,  and  the  mixed-period  parts  of  the  torque 

Now  we  see  that  the  terms  entering  series  (106)  can  be  split  into  three  groups: 

The  terms  with  p  =  h  and  q  =  j  constitute  a  secular  part  of  the  torque,  as  in  such  terms 
vimhj  cancel  with  it*  .  This  M-  and  ^-independent  part  is  furnished  by 


(1) 


2  GM, 


2  R 


2/+1  Z  n  I  oo 

^  i  (/  —  TYl ) !  ^  t  2  ^  t  2 

2-.  (/  +  ,„)!  m  2.  Glpq{e)  kl  Sme'mPq- 

m= 0  '  '  p= 0  q=-o o 

(109) 


1=2 


(2)  The  terms  with  p  =  h  and  q  f=.  j  constitute  a  purely  short-period  part  of  the  torque: 


OO 

t  =  ~yj2 


GM, 


2  R 


2/+1  l  „  ,  l 

2 


1=2 


(/  —  m)\ 

a2,+1  n  (/  +  m)! 

w  ra=0 


ra 


'Y.Flmpd) 


p= o 


oo  oo 

z  z  Gipq(e)  Gipj(e)  k,  sin  [(<7  -  j)  M  -  e,mpq].  (110) 

q=—oo  j=— 00 

(3)  The  terms  with  p  fz  h,  make  a  mixed-period  part  comprised  of  both  short-  and  long- 
period  terms: 

2/+1  1 


_ mixed  -r — ,  2  R  -r — ,  (/  —  m)\  ^ i  ^ — i 

T  m  Z  Z 

t=2  a  m=0  u  T p=o  /j=o 

hjtp 


oo  oo 

z  z  Gihq(e)  Gipj  (e)  kt  sin  [(2/z  -  2/t  +  q  -  j)  M* 

q=—oo  j=—  oo 

+(2h  -  2p)a>*  -  eimpq\. 


(111) 


8.7.3  77;  e  1  =  2  and  1  =  3  terms  in  the  0(i2)  approximation 

For  l  =  2,  index  m  will  take  the  values  0,  1,  2  only.  Although  the  m  =  0  terms  enter  the 
potential,  they  add  nothing  to  the  torque,  because  differentiation  of  (99)  with  respect  to  —9 
furnishes  the  m  multiplier  in  (106).  To  examine  the  remaining  terms,  we  should  consider  the 
inclination  functions  with  subscripts  (Imp)  =  (220),  (210),  (21 1)  only: 

F22od)  =  3  +  0(i2),  F2w(i)  =  -  sin;+  0(i2),  F2u(i)  =  sin i+  O(r),  (112) 

all  the  other  F2mp(i)  being  of  order  0(i 2)  or  higher.  So  for  p  =  h  (i.e.,  both  in  the  secular  and 
short-period  parts)  it  is  sufficient,  in  the  0(i 2)  approximation,  to  keep  only  the  terms  with 
F220 (*') ,  ignoring  those  with  F2l0(i)  and  F2n(i).  We  see  that  in  the  0(i 2)  approximation 
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-  among  the  l  =  2  terms,  both  in  the  secular  and  purely  short-period  parts,  only  the  terms 
with  (Imp)  =  (220)  are  relevant. 

For  p  /:  h,  i.e.,  in  the  mixed-period  part,  the  terms  of  the  leading  order  in  inclination  are: 
Fimp(i)Fimh(i)  =  F2W(i)F2n(i)  and  Fimp(i)F[mh(i)  =  F2\  i  (0  F210(i),  which  happen  to  be 
equal  to  one  another,  and  to  be  of  order  0(i2).  This  way,  in  the  0(i2)  approximation, 

-  the  mixed-period  part  of  the  /  =  2  component  may  be  omitted. 

The  functions  Fhnp  =  F3|0,  F3  l2,  F313,  F32 0,  F32U  F322,  F323,  F33i,  F332 ,  F333  are  of  order 
O  ( i )  or  higher.  The  terms  containing  the  squares  or  cross-products  of  the  these  functions  may 
thus  be  dropped.  Specifically,  the  smallness  of  the  cross-terms  tells  us  that 

-  the  mixed-period  part  of  the  /  =  3  component  may  be  omitted. 

What  remains  is  the  terms  containing  the  squares  of  functions 

^3li(0  =  -^  +  O(i2)  and  F33o(0=  15+O0'2).  (113) 


In  other  words, 

-  among  the  /  =  3  terms,  both  in  the  secular  and  purely  short-period  parts,  only  the  terms 
with  (Imp)  =  (311)  and  (Imp)  =  (330)  are  important  in  the  0(i2)  approximation. 

All  in  all,  for  l  =  2,  3,  the  mixed-period  part  of  the  torque  may  be  dropped  in  the  0(r) 
approximation.  The  secular  and  the  short-period  parts  will  be  developed  up  to  e6,  inclusively. 

8.2  Approximation  for  the  secular  and  short-period  parts  of  the  tidal  torque 

As  we  saw,  the  secular  and  short-period  parts  of  the  torque  may  be  approximated  with 

T  =  T<.2+7,.3+°(e(S/«)9> 

=  T(Ppp)=(220)  +  [T</'»pl=<31!>  +T(/„p)=( 330)]  +  0(€'2,+  0  (f  Wfl)9)’  < 1 14> 

and 

T  =  Tl=1  +Tl=3+0{e(R/a)9) 

=  %,np)=( 220)  +  [^(/mp)=(311)  +  ^(lmp)=( 330)]  +  ° ^  ^  +  0  (6  ( R/° )9)’  (115) 

where  the  /  =  2  and  /  =  3  inputs  are  of  the  order  (R/a)5  and  (R/a)1 ,  accordingly,  while  the 
1  =  4,  5,  . . .  inputs  constitute  0(e  (R/a)9). 

Expressions  for  F (lmp)=(12[S) ,  ^(/mp)=(3n)  3  and  T {lmp)=0W)  arc  Presented  in  the  Elec¬ 
tronic  Supplementary  Material  3  “The  /  =  2  and  1  =  3  terms  of  the  secular  part  of  the 
torque”.  As  an  example,  here  we  provide  one  of  these  components: 
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T 

(lmp)=( 220) 


V  GM2sec  R5  a  6  2^4  e6  ^2  sin  k220_3  |  sgn  (-«  -20) 

{\e2  ~  h  ^  +  7^  g6)  ^  Sin  k220-i 1  Sgn  {n~2d) 

155  A 

—  e  J  k2  sm|62200|sgn  (n  -  6) 

k2  sin  1^2201  I  sgn  On -26) 


■(' 

(: 


1  -5e2  + 


e4  _  155 


49 

T 

'289 

T 

714025 


,  861  4  21975 

e - ~ry  e  +  — — -  e 


’)■ 


/  289  4  1955  A 

(-4  -  e  jk2sm  |62202  I  sgn  (2  n  -  6) 


(116) 


Here  each  term  changes  its  sign  on  crossing  the  appropriate  resonance.  The  change  of  the 
sign  takes  place  smoothly,  as  the  value  of  the  term  goes  through  zero — this  can  be  seen  from 
formula  (91)  and  from  the  fact  that  the  tidal  mode  (Oj  vanishes  in  the  Impq  resonance. 

Expressions  for  T  ,„ni,  T,  and  T.  are  given  in  the  Electronic 

Supplementary  Material  4  “The  l  =  2  and  l  =  3  terms  of  the  short-period  part  of  the  torque”. 
Although  the  average  of  the  short-period  part  of  the  torque  vanishes,  this  part  does  con¬ 
tribute  to  dissipation.  Oscillating  torques  contribute  also  to  variations  of  the  surface  of  the 
tidally-distorted  primary,  the  latter  fact  being  of  importance  in  laser-ranging  experiments. 

The  modes  n  (q  —  j )  of  the  short-period  torque  are  integers  of  n  and  thus  are  commensurate 
with  the  spin  rate  6  near  an  A/S  resonance,  A  and  B  being  integer.  It  may  be  interesting 
to  explore  numerically  the  role  of  this  torque  when  q  —  j  =  1  and  A/B  =  N  is  integer. 
The  possible  role  of  the  short-period  torque  in  the  entrapment  and  libration  dynamics  has 
never  been  discussed  so  far,  as  the  previous  studies  employed  expressions  for  the  tidal  torque, 
which  were  obtained  through  averaging  over  the  period  of  the  secondary’s  orbital  motion. 


9  Marking  the  minefield 

Our  expressions  for  the  secular  and  purely  short-period  parts  of  the  tidal  torque  look  cumber¬ 
some  when  compared  to  the  compact  and  simple  formulae  employed  in  the  literature  hitherto. 
It  will  therefore  be  important  to  explain  why  those  simplifications  are  impractical. 


9. 1  Perils  of  the  conventional  simplification 


Insofar  as  the  quality  factor  is  large  and  the  lag  is  small  (i.e.,  insofar  as  sin  e  can  be  approxi¬ 
mated  with  e),  the  {Imp}  =  {220}  part  of  the  torque  assumes  a  simpler  form: 


(<2>10)_  3  , 

T  =  -  GM1 

1=2  2 


R5a  6k 2  ^ 


G;0q(e)e220q  +  O(e6e)+O(i2e)+O(e\ 
<?=~3  (117) 


where  the  error  0(e3)  originates  from  sine  =  e  +  0(e3). 

The  simplification  conventionally  used  in  the  literature  ignores  the  frequency-dependence 
of  the  Love  number  and  attributes  the  overall  frequency-dependence  to  the  lag.  It  also  ignores 
the  difference  between  the  tidal  lag  e  and  the  lag  in  the  material,  8.  This  way,  the  conventional 
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simplification  makes  the  tidal  lag  e  obey  the  scaling  law  (75b).  At  this  point,  most  authors 
also  set  a  =  —  1 .  Here  we  shall  explore  this  approach,  though  shall  keep  a  arbitrary.  From 
the  formula  1 7 


^ Hmpq 


—Of 

=  £ 


-(<*+!) 

Impq 


(118) 


derived  by  Efroimsky  and  Lainey  (2007)  in  the  said  approach,  we  see  that  the  time  lags  are 
related  to  the  principal-frequency  lag  A r_,_,00  via: 


Qf+1 

(119) 

When  the  despinning  is  still  going  on  and  9  n,  the  corresponding  phase  lags  are: 

o;+l  a 

-  Atlmpq^mpq  =  -At2200Xlmpq  =-*2200(7^)  -  <12°) 

\  ^ Impq  /  \  ^ Impq  / 

which  helps  us  to  cast  the  secular  part  of  the  torque  into  the  following  convenient  form18: 


(<2  >  10)  _ 


1=2 


r  .  /  is  2 

105  4  \  /  i 

( 15  \ 

\-9  1  +  —  e"  + 

—  e4  +  0(e6)]+  n(l  + 

— - 6a  J 

L  V  2 

4  /  V  ' 

^  2  ) 

( 105 

363  \  ,  6  \1 

(t 

— fT“) e  +  °^e 7 

+  0(i-/Q)  +  0(QT 3)  +  0(ae2  Q~1n/9) 

(121a) 


2: 


-9 


{1  +  Ye2)+n  (1  +  (y_6“) g2). 


where  the  overall  factor  reads  as: 


(121b) 


2  = 


^GMseck2 Ar22oo^5  _  3«2MJCC&2 At2200R5 


prim  “{“  Msec)a 

3  n  M^c  k2  R5 


Q220o(Mprim  +  Msec)a  X22OO 


(122) 


17  1 

1 7  sin  €  =  ( ^  )a  >  where  £  is  an  empirical  parameter  of  the  dimensions  of  time,  while  e  is  small  enough, 

so  sin  6  ^  6.  In  combination  with  €,  =  co.  At,  and  y,  =  \co,  I,  this  yields  (118). 

Impq  Impq  impq  Impq  Impq 

18  For  0  2n,  all  the  modes  co are  negative,  so  ^>220q  =  —  ^220 q  '  ^en’  keepin§  in  mind  that  n/0  <$C  1, 

we  process  (120),  for  q  =  1,  as 

a  a 

At  X77nn  ( ^  )  =  -At  2|n  -  B\  (  2ln.~dl  ) 

2200  -200  ^  /2201  J  2200  \  |  -  26  +  3«  \  ) 

=  ~ At22002V-  ”)[1  +  f  7  +  °((^)2) 

L  "  e  J 

and  similarly  for  the  other  values  of  q,  and  then  plug  the  results  into  (117). 
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M prim  and  Msec  being  the  primary’s  and  secondary’s  masses.19  Dividing  (122)  by  the  pri¬ 
mary’s  principal  moment  of  inertia  £  Mprimary  R 2 ,  we  obtain  the  contribution  that  this  com¬ 
ponent  of  the  torque  brings  into  the  angular  deceleration  rate  9 : 


9  =  K.  -  0 


15  ,  105  ,  , 

1  +  y  e2  +  —  e4  +  0(e6) 


+  n 


1  + 


(y -6“)' 


/ 105  363  \ 

I - a  )  i 

V  4  8  ) 


+  _  +  0(c0)  |  ,  +  0tf/Q)  +  o(Q ->)  +  O(ote2 Q~ln/9) 


r  .  /  15  2\ 

(  / 15  ,  \  2\1 

1C 

-9  1  1  +  —  ezJ  +  n 

(l  +  (T-6«)e  ) 

(123a) 

(123b) 


the  factor  K.  being  given  by 
Z 

K.  = 


3  n-  M  k2  At. 


2200 


R 

^  M prim  £  Mprim  (.M prim  3“  Msec)  a 9 

3  n  M2sec  k2 


R  n 


^  Q-t 9qq  Mprim  (Mprim  ^ sec)  t.V'  / 


(124) 


2200 


where  ^  is  a  multiplier  emerging  in  the  expression  £  Mprimary  R2  for  the  primary’s  principal 
moment  of  inertia  (£  =  2/5  for  a  homogeneous  sphere). 

In  the  special  case  of  a  =  —  1 ,  the  above  expressions  enjoy  agreement  with  the  appropriate 
result  stemming  from  the  corrected  MacDonald  model — except  that  our  (121-123)  contain 
A?2'>oO’  Q220O’  ^2200  instead  of  At,  Q,  x  standing  in  formulae  (44-47)  from  Williams  and 
Efroimsky  (2012).  Formula  (123b)  says  that  the  secular  part  of  the  tidal  torque  vanishes  for 


9  —  n  =  —6  n  e 2  a.  (125) 

For  a  =  —  1,  this  coincides  with  the  result  obtained  in  Rodriguez  et  al.  (2008,  eqn.  2.4), 
Correia  et  al.  (2011,  eqn.  20),  and  Williams  and  Efroimsky  (2012,  eqn.  49).  The  coincidence 
however  should  not  be  taken  at  its  face  value,  because  it  is  occasional  or,  possibly  better  to 
say,  exceptional.  Formulae  (121-122)  were  obtained  by  insertion  of  the  expressions  for  the 
eccentricity  functions  and  the  phase  lags  into  (117),  and  by  assuming  that  n  <$;  1 9  \ .  The  latter 
caveat  is  a  crucial  element,  not  to  be  overlooked  by  the  users  of  formulae  (121-122)  and  of 
their  corollary  (123-124)  for  the  tidal  deceleration  rate. 

The  case  of  a  =  —  1  is  special,  in  that  it  permits  derivation  of  (121-125)  without  assum¬ 
ing  that  n  <3C  |0|.  However  for  a  >  —1  the  condition  n  |0|  remains  mandatory,  so 

formulae  (121-124)  become  inapplicable  when  9  reduces  to  values  of  about  several  n. 
Although  formulae  (121a)  and  (123a)  contain  an  absolute  error  0(ae2Q~ln/9),  this  does 

not  mean  that  for  9  close  to  n  the  absolute  error  becomes  O (oie2 Q~l)  and  the  relative  one 

becomes  0(ae 2).  In  reality,  for  9  comparable  to  n,  the  entire  approximation  falls  apart, 
because  formulae  (1 19-120)  were  derived  from  expression  (118),  which  is  valid  for  Q  1 

19  To  arrive  at  the  right-hand  side  of  (122),  we  recalled  that  Ximpq  Atlmpq  =  \elmpq  \  and  that  6^  = 

\elmpq  I  +  0(e3)  =  I €Impq  |  +  0(Q~3),  according  to  formula  (12). 
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only  (unless  a  =  —  1 ).  So  these  formulae  become  inapplicable  in  the  vicinity  of  a  commen- 
surability.  By  ignoring  this  limitation,  one  can  easily  encounter  unphysical  paradoxes.20 

Thence,  in  all  situations,  except  for  the  unrealistic  rheology  a  =  —  1,  limitations  of  the 
approximation  (121-124)  should  be  kept  in  mind.  This  approximation  remains  acceptable 
for  n  <$;  |0|,  but  becomes  misleading  on  approach  to  the  physically-interesting  resonances. 

9.2  An  oversight  in  a  popular  formula 

The  form  in  which  our  approximation  (121-124)  is  cast  may  appear  awkward.  Specifically, 
expression  (124)  for  the  despinning  rate  9  is  written  as  a  function  of  9  and  n,  multiplied  by 
the  factor  1C.  This  form  would  be  reasonable,  were  K.  a  constant.  That  it  is  not — can  be  seen 

from  the  presence  of  the  multiplier  n//2o oo  =  - - - on  the  right-hand  side  of  (124). 

2  |  9  —n\ 

Still,  when  written  in  this  form,  our  result  is  easy  to  juxtapose  with  an  appropriate  for¬ 
mula  from  Correia  and  Laskar  (2004,  2009).  There,  the  expression  for  the  despinning  rate 
looks  similar  to  ours,  up  to  an  important  detail:  the  overall  factor  is  a  constant,  because  it 
lacks  the  said  multiplier  ■  The  multiplier  was  lost  in  those  two  papers,  because  the 
quality  factor  was  introduced  there  as  l/(n  At),  see  the  line  after  formula  (9)  in  Correia  and 
Laskar  (2009).  In  reality,  the  quality  factor  Q  should,  of  course,  be  a  function  of  the  tidal 
frequency  x ,  because  Q  serves  the  purpose  of  describing  the  tidal  damping  at  this  frequency. 
Had  the  quality  factor  been  taken  as  l/(x  At),  it  would  render  the  corrected  MacDonald 
model  (a  =  —  1),  and  the  missing  multiplied  would  be  there.  Being  unphysical,21  the  model 
is  mathematically  convenient,  because  it  enables  one  to  write  down  the  secular  part  of  the 
torque  in  a  compact  form,  avoiding  the  expansion  into  a  series  (Williams  and  Efroimsky 
2012).  The  model  was  pioneered  by  Singer  (1968)  and  employed  by  Mignard  (1979,  1980); 
Hut  (1981)  and  other  authors. 

Interestingly,  in  the  special  case  of  the  3:2  spin-orbit  resonance,  we  have  x  =  n.  Still, 
the  difference  between  x  and  n  in  the  vicinity  of  the  resonance  may  alter  the  probability  of 
entrapment  of  Mercury  into  this  rotation  mode.  The  difference  between  x  and  n  becomes  even 
more  considerable  near  the  other  resonances  of  interest.  So  the  probabilities  of  entrapment 
into  those  resonances  must  be  recalculated. 

10  Conclusions 

The  goal  of  this  paper  was  to  lay  the  ground  for  a  model  of  tidal  entrapment  into  spin-orbital 
resonances.  To  this  end,  we  approached  the  tidal  theory  from  the  first  principles  of  solid-state 
mechanics.  Starting  from  the  expression  for  the  material’s  compliance  in  the  time  domain, 
we  derived  the  frequency-dependence  of  the  Fourier  components  of  the  tidal  torque.  The 
other  torque,  one  caused  by  the  triaxiality  of  the  rotator,  is  not  a  part  of  this  study  and  will 
be  addressed  elsewhere. 

-  We  base  our  work  on  the  Andrade  rheological  model,  providing  arguments  in  favour  of  its 
applicability  to  the  Earth’s  mantle,  and  therefore,  very  likely,  to  other  terrestrial  planets 
and  moons.  The  model  is  also  likely  to  apply  to  the  icy  moons  (Castillo-Rogez  et  al. 


20  For  example,  in  the  case  of  a  >  —  1,  formulae  (118-120)  render  infinite  values  for  A tj  and  on 

crossing  the  commensurability,  i.e.,  when  goes  through  zero. 

21  To  be  exact,  the  model  is  unphysical  everywhere  except  in  the  closest  vicinity  of  the  resonance — see 
formulae  (88-90). 
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2011).  We  have  reformulated  the  model  in  terms  of  a  characteristic  anelastic  timescale 
rA  (the  Andrade  time).  The  ratio  of  the  Andrade  time  to  the  viscoelastic  Maxwell  time, 
£  =  ta/tm,  serves  as  a  dimensionless  free  parameter  of  the  rheological  model. 

The  parameters  rA ,  tm  ,  £  cannot  be  regarded  constant,  though  their  values  may  be  chang¬ 
ing  very  slowly  over  vast  bands  of  frequency.  The  shapes  of  these  frequency-dependencies 
may  depend  on  the  dominating  dissipation  mechanisms  and,  thereby,  on  the  magnitude 
of  the  load,  as  different  damping  mechanisms  get  activated  under  different  loads. 

The  main  question  here  is  whether,  in  the  low-frequency  limit,  anelasticity  becomes  much 
weaker  than  viscosity.  (That  would  imply  an  increase  of  rA  and  f  as  the  tidal  frequency 
X  goes  down.)  The  study  of  ices  under  weak  loads,  with  friction  caused  mainly  by  lat¬ 
tice  diffusion  (Castillo-Rogez  et  al.  2011;  Castillo-Rogez  and  Choukroun  2010)  has  not 
shown  such  a  decline  of  anelasticity.  However,  Karato  and  Spetzler  (1990)  point  out  that  it 
should  be  happening  in  the  Earth’s  mantle,  where  the  loads  are  much  higher  and  damping 
is  caused  mainly  by  unpinning  of  dislocations.  According  to  Karato  and  Spetzler  (1990), 
in  the  Earth,  the  decrease  of  the  role  of  anelasticity  happens  abruptly  as  the  frequency 
falls  below  the  threshold  Xo  ~  1  year- 1 .  We  then  may  expect  a  similar  switch  in  the  other 
terrestrial  planets  and  the  Moon,  though  there  the  threshold  may  be  different  as  it  is  sen¬ 
sitive  to  the  temperature  of  the  mantle.  The  question,  though,  remains  if  this  statement  is 
valid  also  for  the  small  bodies,  in  which  the  tidal  stresses  are  weaker  and  dissipation  may 
be  dominated  by  lattice  diffusion. 

-  We  have  derived  the  frequency  dependencies  of  the  factors  sin  el  emerging  in  the  tidal 
theory.  Naturally,  the  obtained  dependencies  of  these  factors  upon  the  tidal  frequency 
Xjmpq  (°r>  more  exactly,  upon  the  tidal  mode  u>im )  mimic  the  frequency-dependence  of 
the  imaginary  part  of  the  complex  compliance.  They  scale  as  ~  x  “  with  0  <  a  <  1,  at 
higher  frequencies;  and  as  ~  x  1  at  lower  frequencies.  However  in  the  zero-frequency 
limit  the  factors  kt  sin  e:  demonstrate  a  behaviour  inherent  in  the  tidal  lagging  and  absent 
in  the  case  of  lagging  in  a  sample:  in  a  close  vicinity  of  the  zero  frequency,  these  factors 
(and  the  appropriate  components  of  the  tidal  torque)  become  linear  in  the  frequency.  This 
way,  sin  e;  first  reaches  a  finite  maximum,  then  decreases  continuously  to  nil  as  the 
frequency  approaches  to  zero,  and  then  changes  its  sign.  So  the  resonances  are  crossed 
continuously,  with  neither  the  tidal  torque  nor  the  tidal  force  diverging  there.  For  example, 
the  leading  term  of  the  torque  vanishes  at  the  synchronous  orbit. 

The  continuous  traversing  of  resonances  was  pointed  out  in  a  heuristic  manner  by  Efroim¬ 
sky  and  Williams  (2009).  Now  we  have  derived  this  result  directly  from  the  expression 
for  the  compliance  of  the  material  of  the  rotating  body.  Our  treatment  however  has  a  flaw: 
the  frequency,  below  which  the  factors  k:  sin  e,  and  the  appropriate  components  of  the 
torque  change  their  frequency-dependence  to  linear,  is  very  low  (lower  than  1 0  10  Hz,  if 
we  take  our  formulae  literally).  The  reason  for  this  is  that  we  kept  using  the  known  value 
of  the  Maxwell  time  xM  all  the  way  down  to  the  zero  frequency.  Possible  changes  of  the 
Maxwell  time  in  the  zero-frequency  limit  may  broaden  the  region  of  linear  dependence. 

-  We  have  offered  an  explanation  of  the  “improper”  frequency-dependence  of  dissipation  in 
the  Moon,  discovered  by  LLR.  The  main  point  of  our  explanation  is  that  the  LLR  measures 
the  tidal  dissipation  whose  frequency-dependence  is  different  from  that  of  the  seismic 
dissipation.  Specifically,  the  “wrong”  sign  of  the  exponential  in  the  power  dissipation 
law  may  indicate  that  the  frequencies  at  which  tidal  friction  was  observed  were  below  the 
frequency  at  which  the  lunar  kj  sin  62  has  its  peak.  Taken  the  relatively  high  frequencies 
of  observation  (corresponding  to  periods  of  order  month  to  year),  this  explanation  can  be 
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accepted  only  if  the  lunar  mantle  has  a  low  mean  viscosity.  This  may  be  the  case,  taken 
the  presumably  high  concentration  of  the  partial  melt  in  the  low  mantle. 

-  We  have  developed  a  detailed  formalism  for  the  tidal  torque,  and  have  singled  out  its 
oscillating  component. 

The  studies  of  entrapment  into  spin-orbit  resonances,  performed  in  the  past,  took  into 
account  neither  the  afore-mentioned  complicated  frequency-dependence  of  the  torque  in 
the  vicinity  of  a  resonance,  nor  the  oscillating  part  of  the  torque. 
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Appendix 

A  Interconnection  between  the  quality  factor  and  the  phase  lag 


The  power  P  exerted  by  a  tide-raising  secondary  on  its  primary  can  be  written  as 

P  =  -  J  p  V  •  VW  d3x  (126) 

p,  V,  and  W  signifying  the  density,  velocity,  and  tidal  potential  in  the  small  volume  d3 x  of 
the  primary.  Since  V  •  (p  V)  +  ^  =  0,  the  dot  product  takes  the  form  of 

p  v.  vw  =  v  •  (p  v  w)  -p  w  v  •  v-  v  w  vp.  (127) 

Under  the  realistic  assumption  of  the  primary’s  incompressibility,  the  term  with  V  •  V  may 
be  omitted.  To  get  rid  of  the  term  with  Vp,  one  has  to  accept  a  much  stronger  approximation 
of  the  primary  being  homogeneous.  Then  the  power  will  be  rendered  by 

P  =  -  V  ■  (p  V  W)d3x  =  -  pW\-ndS,  (128) 


n  being  the  outward  normal  and  dS  being  an  element  of  the  surface  area  of  the  primary.  This 
expression  for  the  power  (pioneered,  probably,  by  Goldreich  1963)  enables  one  to  calculate 
the  work  through  radial  displacements  only,  in  neglect  of  horizontal  motion. 

Denoting  the  radial  elevation  with  f ,  we  write  the  power  per  unit  mass,  V  =  P/M,  as: 

A  harmonic  external  potential 

W  =  Wo  cos {colmpq  t),  (130) 

applied  at  a  point  of  the  primary’s  surface,  will  elevate  this  point  by 

?  =  /;2  ~r  cos  [a)lmpq  t  -  elmpqj  =  h2  — ^  cos  [wlmpq  t  -  wlmpq  A tlmpqJ,  (131) 

D  O 

with  g  being  the  surface  gravity  acceleration,  and  h2  denoting  the  Love  number. 
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In  formula  (131),  ojlmpq  is  one  of  the  modes  (101)  showing  up  in  the  Darwin-Kaula 
expansion  (99).  The  quantity  elmpq  =  a>lmpc/  A tlm  is  the  corresponding  phase  lag,  while 
A is  the  positively  defined  time  lag  at  this  mode.  Although  the  tidal  modes  (o{  can 
assume  any  sign,  both  the  potential  W  and  elevation  t,  can  be  expressed  via  the  positively 
defined  forcing  frequency  Ximpq  =  I  wimpq  I  and  the  absolute  value  of  the  phase  lag: 


W  =  Wq  cos(x  t), 

W 

S  =  h2—  cos(x  t  -  |  6  |), 


(132) 

(133) 


subscripts  Impq  being  dropped  here  and  hereafter  for  brevity. 

The  vertical  velocity  of  the  considered  element  of  the  primary’s  surface  will  be 

d ?  ,  W„  .  ,  .  ..  W0  .  .  .  . 

-r  =  ~h2X —  sin(xf  -  |e|)  =  -li2X —  (sin  xf  coski  -  cosxt  sin|e|) . 

dt  g  g 

Introducing  the  notation  A  =  h2  ,  we  write  the  power  per  unit  mass  as 

(- 


,  dW\  dt 

V  =  (  — —  \  —  =  Ax  cos(x  t)  sin (xt  -  |e|), 


(134) 


(135) 


and  express  the  work  w  per  unit  mass,  performed  over  a  time  interval  (to,  t ),  as: 
t  xt 

cos(x  t)  sin(x?  -  |e|)r/(x  0 


w\t  =  V  dt  =  A 


to  Xto 

xt  xt 

=  Acos|e|  /  cosz  sin  zdz  —  A  sin  |e|  j  cos 2  zdz 

x>o 


X‘ 

x'o 

A  t 

=  -  ----  [cos(2xt  -  |e|)  +2  x  t  sink 


to 


(136) 


Being  cyclic,  the  first  term  in  (136)  renders  the  elastic  energy  stored  in  the  body.  The  second 
term,  being  linear  in  time,  furnishes  the  energy  damped.  This  clear  interpretation  of  the  two 
terms  was  offered  by  Stan  Peale  (201 1,  personal  communication). 

The  work  over  a  time  period  T  =  2n/x  is  equal  to  the  energy  dissipated  over  the  period: 


t=T 

w|,=0  =  AEcycle  =  ~ A  7T  sin|e|. 


(137) 


The  peak  work  is  obtained  over  the  time  span  from  k  to  |e|  and  assumes  the  value 

(work) 
peak 


A  r 

(n  ,  ,M 

—  cos  k  —  sin  k 

2  L 

V  2  7J 

(138) 


whence  the  appropriate  quality  factor  is  given  by: 

—A  E 


Q-1  ,  = 

work 


‘‘cycle 


tan  lei 


2  j r  E 


(work) 

peak 


1  -  (f  -  kl) 


tan  kl 


(139) 


To  calculate  the  peak  energy  stored  in  the  body,  we  would  note  that  the  first  term  in  (136) 
is  maximal  when  taken  over  the  span  from  x  t  =  ;r/ 4  +  \e\/2  through  x  /  =  3^/4  +  kl/2: 


_( energy ) 
^ peak 


A 

2’ 


(140) 
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and  the  corresponding  quality  factor  is: 


Q~l 

energy 


Cyrfe 


2  jt  E 


( energy ) 
peak 


=  sin  |e|. 


(141) 


Goldreich  (1963)  suggested  to  employ  the  span  x  t  =  (0,  tt/4).  The  absolute  value  of 
the  resulting  power,  denoted  by  Goldreich  as  E* ,  is  equal  to 


£*=y  cos|e|  (142) 

and  is  neither  the  peak  value  of  the  energy  stored,  nor  that  of  the  work  performed.  Goldreich 
(1963)  however  employed  it  to  define  a  quality  factor,  which  we  shall  term  Q(}0idreich  • 


Q 


-i 

Goldreich 


-A  E 


cycle 


2tt  E* 


=  tan  |e|. 


(143) 


The  quality  factor  Qenergy  defined  through  (141)  is  preferable,  as  the  expansion  of  tides 
contains  terms  proportional  to  ki(ximpq)  sin€;(//m;?„)-  Since  the  long-established  tradition 
suggests  to  substitute  sin  |e|  with  1  / Q,  it  is  advisable  to  define  the  Q  exactly  as  (141),  and 
also  to  call  it  Qi,  to  distinguish  it  from  the  seismic  quality  factor  (Efroimsky  2012). 


B  The  correspondence  principle  (elastic-viscoelastic  analogy) 


B.l  The  correspondence  principle,  for  nonrotating  bodies 


While  the  static  Love  numbers  depend  on  the  static  rigidity  p  through  (3),  it  is  not  imme¬ 
diately  clear  if  a  similar  formula  interconnects  also  ki(x)  with  p(x)-  To  understand  why 
and  when  the  relation  should  hold,  recall  that  formulae  (3)  originate  from  the  solution  of  a 
boundary- value  problem  for  a  system  incorporating  two  equations: 


aBv  —  2^  ut 


0  = 


da, 


flv 


dx„ 


dp 

dxo 


-  P 


3  (W  +  U) 
dxK 


(144a) 


the  latter  being  simply  the  equation  of  equilibrium  written  for  a  static  viscoelastic  medium, 
in  neglect  of  compressibility  and  heat  conductivity.  The  notations  a^v  and  u^v  stand  for  the 

deviatoric  stress  and  strain,  p  =  —  ^  Sp  S  is  the  pressure  (set  to  be  nil  in  incompressible 
media),  while  W  and  U  are  the  perturbing  and  perturbed  potentials.  By  solving  the  system, 
one  arrives  at  the  static  relation  (7/  =  k/  W/,  with  the  customary  static  Love  numbers  ki 
expressed  via  p,  R ,  and  p.  by  (3). 

Generalise  (144a-144a)  to  the  case  of  time-dependent  loading  of  a  nonrotating  body: 


S  =  2£U,  (145a) 

pu  =  VS—  Vp  — pV(W  +  (7)  (145b) 

or,  in  terms  of  components: 

apv  =  2A  upv , 

d°fSv  dp  3  (W  +  U) 

“  dxv  dXp  dXp 
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In  the  frequency  domain,  this  will  look: 


V(x) 

PX2«/jv,(X) 


2  ft(x)  (x)> 

9y(x)  3p(X)  9[W(/)  +  t/(x)] 

3jcu  3x^ 


(147a) 

(147b) 


where  a  bar  denotes  a  spectral  component  for  all  functions  except  p. — recall  that  [i  is  a 
spectral  component  not  of  the  kernel  ^(r)  but  of  its  time-derivative 

Unless  the  frequencies  are  extremely  high,  we  can  neglect  the  body-fixed  acceleration  term 
X2  in  the  second  equation,  in  which  case  our  system  of  equations  for  the  spectral 

components  will  mimic  (144).  Thus  we  arrive  at  the  so-called  correspondence  principle  (also 
known  as  the  elastic-viscoelastic  analogy),  which  maps  a  solution  of  a  linear  viscoelastic 
boundary-value  problem  to  a  solution  of  a  corresponding  elastic  problem  with  the  same  initial 
and  boundary  conditions.  As  a  result,  the  algebraic  equations  for  the  Fourier  (or  Laplace) 
components  of  the  strain  and  stress  in  the  viscoelastic  case  mimic  the  equations  connecting 
the  strain  and  stress  in  the  appropriate  elastic  problem.  So  the  viscoelastic  operational  moduli 
/x(x)  or  J(x)  obey  the  same  algebraic  relations  as  the  elastic  parameters  /x  or  J . 

In  the  literature,  there  is  no  consensus  on  the  authorship  of  this  principle.  For  example, 
Haddad  (1995)  mistakenly  attributes  it  to  several  authors  who  published  in  the  1950s  and 
1960s.  In  reality,  the  principle  was  pioneered  almost  a  century  earlier  by  Darwin  (1879),  for 
isotropic  incompressible  media.  The  principle  was  extended  to  more  general  types  of  media 
by  Biot  (1954,  1958),  who  also  pointed  out  some  limitations  on  its  applicability. 


B.2  The  correspondence  principle,  for  rotating  bodies 


Consider  a  body  of  mass  Mprjm ,  which  is  spinning  at  a  rate  (o  and  is  also  performing  some 
orbital  motion  (for  example,  is  orbiting,  with  its  partner  of  mass  Msec,  around  their  mutual 
centre  of  mass).  Relative  to  some  inertial  coordinate  system,  the  centre  of  mass  of  the  body 
is  located  at  xCM ,  while  a  small  parcel  of  its  material  is  positioned  at  x.  Relative  to  the  centre 
of  mass  of  the  body,  the  parcel  is  located  at  r  =  x  —  xCM .  The  body  being  deformable,  we 
can  decompose  r  into  its  average  value,  ro,  and  an  instantaneous  displacement  u: 


X  =  XCM+r 

r  =  r0  +  u 


■  x  =  xrM  +  r0  +  u. 


(148) 


Denote  with  D/Dt  the  time-derivative  in  the  inertial  frame.  The  symbol  d/dt  and  overdot 
will  be  reserved  for  the  time-derivative  in  the  body  frame,  so  dro/dt  =  0.  Then 

Dr  dr  D2  r  d2r  dr 

—  = - E  to  x  r  and  — =  — =-  +  2co  x - b  co  x  (at  x  r)  +  co  x  r.  ( 149) 

Dt  dt  Dt2  dt2  dt 

Together,  the  above  formulae  result  in 


D2x 

H)!2 


D2x  D2 r 

^  acm  ,  _  1 

Dt 2  Dt 2 


Dzxr 


Dt2 


d2r  dr 

-\ - =-  +2(0  x - b  «  x  (<nxr)  +  toxr 

dtA  dt 


D2xn 


Dt 2 


d2  u  du 

-| - =-  +  2oo  x - b  (o  x  (a)  x  r)  +  (o  x  r. 

dt 2  dt 


(150) 


The  equation  of  motion  for  a  small  parcel  of  the  body’s  material  will  read  as 


D-x 

P  2  =  —  V p  +  F self  +  F ext, 


(151) 
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where  FfA-r  is  the  exterior  gravity  force  per  unit  volume,  while  F self  is  the  “interior”  gravity 
force  per  unit  volume,  i.e.,  the  self-force  wherewith  the  rest  of  the  body  is  acting  upon  the 
selected  parcel  of  medium.  Insertion  of  (150)  in  (151)  furnishes: 


nzx 

aca 
Dt 2 


+  U  +  2a>  x  u  +«  x(«xr)  +  «xr 


—  VS  —  V  p  +  F  se]f  +  Fe 


(152) 


At  the  same  time,  for  the  primary  body  as  a  whole,  we  can  write: 


D2xr 


Dt 2 


(153) 


the  integration  being  carried  out  over  the  volume  V  of  the  primary.  (Recall  that  Fext  is  a  force 
per  unit  volume.)  Combined  together,  the  above  two  equations  will  result  in 

pj\j  +  2o>x  u+wx  (<Dxr)  +  wxrJ 

=  VS  -  V p  +  Fself  +  Fext  -  [  Fext  d\.  (154) 

Mprim  J 
v 


For  a  spherically-symmetrical  (not  necessarily  radially-homogeneous)  body,  the  integral  on 
the  right-hand  side  removes  the  Newtonian  part  of  the  force,  leaving  the  harmonics  intact: 


F„,  -  — —  /  Fext  d3r  =  p  ]T  Vff,,  (155) 

M prim  J 
y 

where  the  harmonics  are  given  by 

.  GM'er  /  r  \! 

Wt(  r,  r  )  = - —  (—  j  Pi  (cos  y),  (156) 

r*  being  the  vector  pointing  from  the  centre  of  mass  of  the  primary  to  that  of  the  secondary, 
and  y  being  the  angle  between  r  and  r* ,  subtended  at  the  centre  of  mass  of  the  primary. 

In  reality,  a  tiny  extra  force  T,  the  tidal  force  per  unit  volume,  is  left  over  due  to  the  body 
being  slightly  distorted: 


Fgjrr 


M 


prim 


oo 

1=2 


(157) 


where  T  is  the  density  multiplied  by  the  average  tidal  acceleration  experienced  by  the  body 
as  a  whole.  In  neglect  of  T ,  we  arrive  at 


OO 

p[u+  2o>x  u  +  a)  x  («  x  r)  +  U  x  r j  =  VS  —  Vp  —  p  ^V(f//  +  lV/).  (158) 

1=2 

Here,  to  each  disturbing  term  of  the  exterior  potential,  W/,  corresponds  a  term  U\  of  the 
self-potential,  the  self-force  thus  being  expanded  into  F seif  =  —  Vt//. 

Equation  (158)  could  as  well  have  been  derived  in  the  body  frame. 

Denoting  the  tidal  frequency  with  x ,  we  see  that  the  terms  on  the  left-hand  side  have  the 
order  of  p  /2  u,  p  to  /  u,  p  co2  r ,  and  p  /  cor,  correspondingly.  In  realistic  situations,  the  first 
two  terms,  thus,  can  be  neglected,  and  we  end  up  with 
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oo 

0  =  VS  —  Vp  —  p  ~y  V  {Ui  +  Wi)  —  poo  x  (w  x  r)  —  p  w  x  r,  (159) 

1=2 

the  term  —  V  p  vanishing  in  an  incompressible  media. 


B.3  The  centripetal  term  and  the  zero-degree  Love  number 


The  centripetal  term  in  (159)  can  be  split  into  a  purely  radial  part  and  a  part  that  can  be 
incorporated  into  the  W2  term  of  the  tide-raising  potential,  as  was  suggested  by  Love  (1909, 
1911).  Introducing  the  colatitude  (p'  through  cos (p '  =  •  -p,  we  can  write: 


to  x  (w  x  r)  =  co  (to  ■  r)  —  r  w2  =  V 


2  „2 


-  (to  •  r)  —  to  r 
2  2 


=  V 


-  to"  r"  (cos 2  4>'  ~  1) 


The  definition  P2(coscp')  =  A  (3  cos 2  (p1  —  l)  renders:  cos2  cp'  =  %  P2(coscp')  +  A 


whence: 


to  x  (to  x  r)  =  V 


2  2 

-  w  r“ 


[P2  (coscp')  -  l] 


3’ 
(160) 


We  see  that  the  centripetal  force  splits  into  a  second-harmonic  and  purely-radial  parts: 

—  pw  x  (to  x  r)  =  —  V  |j: ,  to2r2  P2  (cos0')j  +  V  |^/i-w2r2j,  (161) 

where  we  assume  the  body  to  be  homogeneous.  The  second-harmonic  part  can  be  incorpo¬ 
rated  into  the  external  potential.  The  response  to  this  part  will  be  proportional  to  the  degree-2 
Love  number  k2. 

The  purely  radial  part  of  the  centripetal  potential  generates  a  radial  deformation.  This 
part  of  the  potential  is  often  ignored,  the  associated  deformation  being  tacitly  included 
into  the  equilibrium  shape  of  the  body.  Compared  to  the  main  terms  of  the  equation  of 
motion,  this  radial  term  is  of  the  order  of  1 0  '  for  the  Earth,  and  is  smaller  for  most  other 
bodies.  As  the  rotation  variations  of  the  Earth  are  of  the  order  of  10~5,  this  term  leads 
to  a  tiny  change  in  the  geopotential  and  to  an  associated  displacement  of  the  order  of  a 

•  9? 

micrometer. 

However,  for  other  rotators  the  situation  may  be  different.  For  example,  in  Phobos,  whose 
libration  magnitude  is  large  (about  1  degree),  the  radial  term  may  cause  an  equipotential- 
surface  variation  of  about  10  cm.  This  magnitude  is  large  enough  to  be  observed  by  future 
missions  and  should  be  studied  in  more  detail. 23  The  emergence  of  the  purely  radial  defor¬ 
mation  gives  birth  to  the  zero-degree  Love  number  (Dahlen  1976;  Matsuyama  and  Bills 
2010).  Using  Dahlen’s  results,  Yoder  (1982,  eqns  21-22)  demonstrated  that  the  contribution 
of  the  radial  part  of  the  centripetal  potential  to  the  change  in  mean  motion  of  Phobos  is  about 
3%,  which  is  smaller  than  the  uncertainty  in  our  knowledge  of  Phobos’  k2/ Q.  It  should  be 
mentioned,  however,  that  the  calculations  by  Dahlen  (1976)  and  Matsuyama  and  Bills  (2010) 
were  performed  for  steady  (or  slowly  changing)  rotation,  and  not  for  libration.  So  Yoder’s 
application  of  Dahlen’s  result  to  Phobos  requires  extra  justification. 

What  is  important  for  us  here  is  that  the  radial  term  does  not  interfere  with  the  calcula¬ 
tion  of  the  Love  numbers  of  non-zero  degree.  Being  independent  of  the  longitude,  this  term 
generates  no  tidal  torque  either,  provided  the  obliquity  is  neglected. 


22  Tim  Van  Hoolst,  private  communication. 
-  '  Tim  Van  Hoolst,  private  communication. 
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B. 4  The  toroidal  term 

The  inertial  term  — p  wxrin  the  equation  of  motion  (159)  can  be  cast  into  the  form 

—  p  d)  x  r  =  p  r  x  V(<h  ■  r),  (162) 

whence  we  see  that  this  term  is  of  a  toroidal  type.  Being  almost  nil  for  a  despinning  primary, 
this  force  becomes  important  for  a  librating  object. 

In  spherically-symmetric  bodies,  the  toroidal  force  (162)  generates  toroidal  deformation 
only.  The  latter  produces  neither  radial  uplifts  nor  variations  of  the  gravitational  potential. 
Hence  its  presence  does  not  influence  the  expressions  for  the  Love  numbers  associated  with 
vertical  displacement  )  or  the  potential  (kt ).  As  this  deformation  yields  no  change  in  the 
gravitational  potential  of  the  tidally-perturbed  body,  it  generates  no  tidal  torque  either.  Being 
divergence-free,  this  deformation  entails  no  contraction  or  expansion,  i.e.,  it  is  purely  shear. 
Still,  this  deformation  contributes  to  dissipation.  Besides,  since  the  toroidal  forcing  yields 
toroidal  deformation,  it  can,  in  principle,  be  associated  with  a  “toroidal  Love  number”. 

To  estimate  the  dissipation  caused  by  the  toroidal  force,  Yoder  (1982)  introduced  an  equiv¬ 
alent  effective  torque.  He  pointed  out  that  the  force  gets  important  when  the  magnitude  of  the 
physical  libration  is  comparable  to  that  of  the  optical  libration.  According  to  Yoder  (1982), 
the  toroidal  force  contributes  to  the  change  of  Phobos’  mean  motion  about  1.6%,  which  is 
less  than  the  input  from  the  purely  radial  part. 

C  The  Andrade  and  Maxwell  models  at  different  frequencies 

C. l  Response  of  a  sample  obeying  the  Andrade  model 

Within  the  Andrade  model,  the  tangent  of  the  phase  lag  demonstrates  the  so-called  “elbow 
dependence”.  At  high  frequencies,  the  tangent  of  the  lag  obeys  a  power  law  with  an  exponent 
equal  to  —a,  where  0  <  a  <  1.  At  low  frequencies,  the  tangent  of  the  lag  once  again  obeys  a 
power  law,  this  time  though  with  an  exponent  —(1  —  a).  This  model  fits  well  the  behaviour 
of  ices,  metals,  silicate  rocks,  and  many  other  materials. 

However  the  applicability  of  the  Andrade  law  may  depend  on  the  intensity  of  the  load 
and,  accordingly,  on  the  damping  mechanisms  involved.  Situations  are  known,  when,  at  low 
frequencies,  anelasticity  becomes  much  less  efficient  than  viscosity,  and  the  Andrade  model 
approaches  the  Maxwell  model. 

C.1.1.  The  high-frequency  band 

At  high  frequencies,  expression  (85b)  gets  simplified.  In  the  numerator,  the  term  with  z~a 

l 

dominates:  z~a  ~S>  z  1  f ,  which  is  equivalent  to  z  f  1  a  .  In  the  denominator,  the  constant 

l 

t  #  j _ ^ 

term  dominates:  1  z~a,  or  simply:  z  1.  To  know  which  of  the  two  conditions,  z  f 

or  z  1 ,  is  stronger,  we  recall  that  at  high  frequencies  anelasticity  beats  viscosity.  So  the 
a-term  in  (81)  is  large  enough.  In  other  words,  the  Andrade  timescale  rA  should  be  smaller 
(or,  at  least,  not  much  higher)  than  the  viscoelastic  time  rM .  Accordingly,  at  high  frequen¬ 
cies,  f  is  smaller  (or,  at  least,  not  much  higher)  than  unity.  Hence,  within  the  high-frequency 

l 

1  — OL 

band,  either  the  condition  z  1  is  stronger  than  z  K  or  the  two  conditions  are  about 
equivalent.  This,  along  with  (86)  and  (87)  enables  us  to  write: 
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tan<5«(xrA)  sin  V(a  +  1)  for  x  »  rAl=zM1?  '•  (163) 

The  tangent  being  small,  the  expression  for  sin  S  looks  identical: 

sin<5«(xrA)  sin(^)r(a+l)  for  x  »  r"1  =  r"1^1.  (164) 


C.1.2.  The  intermediate  region 

In  the  intermediate  region,  the  behaviour  of  the  phase  lag  S  depends  upon  the  frequency- 

dependence  of  f .  For  example,  if  there  happens  to  exist  an  interval  of  frequencies  over  which 
l 

^  t  J  _ Ql  ^  f 

the  conditions  1  z  f  are  obeyed,  then  over  this  interval  we  shall  have:  1  z~a 
and  z~“  f  z^1 .  Applying  these  inequalities  to  (85b),  we  see  that  over  such  an  interval  of 
frequencies  tan  S  will  behave  as  z~2“  tan  I  ^ -  j . 


C.1.3  The  low-frequency  band 


At  low  frequencies,  the  term  z  1  f  becomes  leading  in  the  numerator  of  (85b):  z  01  <SC  z  1  f , 
l 

i  i  2 _ (x  ,  . 

which  requires  z  f  .  In  the  denominator,  the  term  with  z~“  becomes  the  largest: 
1  -C  z~“,  whence  z  <SC  1.  Since  at  low  frequencies  the  viscous  term  in  (81)  is  larger 
than  the  anelastic  term,  we  expect  that  for  these  frequencies  f  is  larger  (at  least,  not  much 

smaller)  than  unity.  Thence  the  condition  z  -C  1  becomes  sufficient.  Its  fulfilment  ensures 
l 

1  — Cl 

the  fulfilment  of  z  <g  (  .  Thus  we  state: 


tan ,5  »  (x  rA) 


-(!-<*) 


•{*£)  r(«  +  i) 


for  x  «  \'=rMlS  l- 


(165) 


The  appropriate  expression  for  sin  S  will  be: 


sin 5  «  1  -  O  ((x  rA)2<1-“)  ?"2)  for  X  «  rf1  =  t r*.  (166) 

It  would  be  important  to  emphasise  that  the  threshold  zfl  =  zM  1  standing  in  (163) 
and  (164)  is  different  from  the  threshold  z~l  =  z^1 1,  showing  up  in  (165)  and  (166), 
even  though  these  two  thresholds  are  given  by  the  same  expression.  The  reason  for  this  is 
that  the  timescales  rA  and  zM  are  not  fixed  constants.  While  the  Maxwell  time  is  likely  to  be 
a  very  slow  function  of  the  frequency,  the  Andrade  time  may  undergo  a  faster  change  over 
the  transitional  region:  rA  must  be  larger  than  xM  at  low  frequencies  (so  anelasticity  yields 
to  viscosity),  and  must  become  shorter  than  or  of  the  order  of  tm  at  high  frequencies  (so 
anelasticity  becomes  stronger).  This  way,  the  threshold  zf1  standing  in  (165-166)  is  lower 
than  the  threshold  zfl  standing  in  (163-164).  The  gap  between  these  thresholds  is  the  region 
intermediate  between  the  two  pronounced  power  laws  (163)  and  (165). 


C.1.4  The  low-frequency  band:  a  special  case,  the  Maxwell  model 

Suppose  that,  below  some  threshold  Xo ,  anelasticity  quickly  becomes  much  less  efficient  than 
viscosity.  This  would  imply  a  steep  increase  of  f  (equivalently,  of  rA)  at  low  frequencies. 
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Then,  in  (85b),  we  shall  have:  1  »  z  “  and  z  “  <SC  fz  *.  This  means  that,  for  frequencies 
below  Xo>  the  tangent  will  behave  as 

tanS  ss  z_1  C  =  (x  for  x  «  (16V) 

the  well-known  viscous  scaling  law  for  the  lag. 

The  study  of  ices  and  minerals  under  weak  loads  (Castillo-Rogez  et  al.  2011;  Castillo- 
Rogez  and  Choukroun  20 1 0)  has  not  shown  such  an  abrupt  vanishing  of  anelasticity.  However, 
Karato  and  Spetzler  (1990)  point  out  that  this  should  be  happening  in  the  Earth’s  mantle, 
where  the  loads  are  much  higher  and  anelasticity  is  caused  by  unpinning  of  dislocations. 


C. 2  The  behaviour  of  (x)  |  sine;(x)  =  —  I/n[/q(x)]  within  the  Andrade  and  Maxwell 

models 

As  we  explained  in  Sect.  4.1,  products  (Ximpq  )  sin  (Ximpq )  enter  the  Impq  term  of  the 
Darwin-Kaula  series  for  the  tidal  potential,  force,  and  torque.  Hence  the  importance  to  know 
the  behaviour  of  these  products  as  functions  of  the  tidal  frequency  Ximpq  ■ 

C.2.1  Prefatory  algebra 


It  ensues  from  (64)  that 


k,(x)  = 


3  {Ke[J(X)])2  +  {lm[J(x)])2+A,  JUe[J(x)]  +i  A,  J  lm[J(x)\ 


2(1-1) 


(TZe[J(x)]  +  A,  j)2  +  (lm[J(X)])2 


(168) 


whence 

\k,(x)\  sin  ej  (x ) 


-Xm  [fc,(x)] 


_3 _ -A,Jlm[J(x)] _ 

2  (/  “  1}  {Tie  [  J(x)]  +  A,  J)2  +  (Xm  [/(x)])2  ’ 

(169) 


7  =  7(0)  =  l/p  =  l//r(0)  being  the  unrelaxed  compliance  (the  inverse  of  the  unrelaxed 
shear  modulus  p).  For  an  Andrade  material,  the  compliance  7  in  the  frequency  domain  is 
rendered  by  (82).  Its  imaginary  and  real  parts  are  given  by  (83-84).  It  is  then  easier  to  rewrite 
(169)  as 

3  A, 

\k,(x)\  sme,(x)  =  2(/  _  jy 

?  z^1  +z_“  sin(^)  I>+  1) 

^  2  2  ’ 

[a,  +  1  +  z~“  cos  r(a  +  1)]“  +  [?  z~l  +  z"“  sin  T(1  +  a)]" 

(170) 


where 


z  =  X  rA  =  x  rM  $ 


(171) 


and 


(172) 
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For /l  ->  0,  i.e.,forrA  ->  oo,  (170)  coincides  with  the  appropriate  expression  for  a  spherical 
Maxwell  body. 

C.2.2  The  high-frequency  band 

Within  the  upper  band,  the  term  with  z~“  dominates  the  numerator,  while  At  +  1  dominates 

l 

the  denominator.  The  dominance  of  z~“  in  the  numerator  requires  that  z  t;  ,  which  is: 

a 

.  I— a 

XtM»? 

—  \/a 

The  dominance  of  A?  +  1  in  the  denominator  means:  z  (At  +  1)  ,  which  is  the  same 

as 

i  —l/a 

(ai  +  1  > 

The  dominance  of  At  +  1  also  demands  that  f  z^1  <SC  (A:  +  1),  which  is: 

X  »  (Ai  +  1)  *■ 

Simultaneous  fulfilment  of  all  these  conditions  entails: 

lfy(x)l  sin€,(/)^2(/3-  i)  (A  +1)2  sin(~r)F(0;+  l^~a{T«Xya,  for  X  »  Xhi  > 

(173a) 

where  the  boundary  between  the  high  and  intermediate  frequencies,  xHi  >  *s  defined  as  the 

OL 

i  1 — a  |  i  \  /  a  *  -i 

quantity  largest  among  x ~  t;  ,  r ~  f  (Aj  +  1)  ,  and  xM  (A,  +  1)  . 

For  small  terrestrial  objects,  At  )$>  1  (Efroimsky  2012).  As  at  high  frequencies  anelas- 
ticity  beats  viscosity,  f  may  be  of  order  unity  or  slightly  less  than  unity  (not  by  orders  of 

a 

magnitude,  likely).  Hence  for  small  bodies  and  small  terrestrial  planets  xHI  =  t“*  f  a  ■ 
For  solid  bodies  much  larger  than  the  Earth,  A;  <$;  1  (see  Efroimsky  2012).  Hence  for 
large  terrestrial  exoplanets  Xhi  =  f-1. 

C.2.3  The  intermediate  band 

Within  the  intermediate  band,  the  term  f  z^1  takes  over  in  the  numerator,  while  A;  +  1 

still  dominates  in  the  denominator.  The  dominance  of  f  z_1  in  the  numerator  implies  that 
l 

J _ #  i  f 

z  <3C  f  ,  which  is  equivalent  to 

a 

.  1  —O' 

X  rM  «  £ 

The  dominance  of  At  +  1  in  the  denominator  requires,  as  we  just  saw: 

-i  ,  -l/a 

X  tm  ^  f  (At  +  1) 

and 

X  r*f  »  (A;  +  1)  1  - 
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We  are  considering  the  band  where  the  efficiency  of  viscosity  is  either  comparable  or  exceeds 
that  of  anelasticity.  Accordingly,  f  is  about  or  larger  than  unity.  Then,  of  the  two  above-written 
inequalities,  the  latter  is  stronger,  for  any  A; .  Thus  we  obtain: 


a 


3  A, 


(tmX)  \  for  TMlt;  1  “  »X»rM1(A/  +  l)  *. 

(173b) 


lfy(x)l  sine/(x) 


2  (/  —  1)  (A,  +  1)2 


C.2.4  The  low-frequency  band 

For  frequencies  lower  than  (A/  +  l)-1,  a  self-gravitating  Andrade  body  renders  the  same 
frequency-dependency  as  a  self-gravitating  Maxwell  body: 
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